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FOREWORD 

This  report  describes  a  contractual  work  effort  conducted 
for  the  General  Electric  Company,  Aircraft  Engine  Group  under 
Purchase  Order  Mo.  200-FBA-17K-47844  which  is  a  subcontract 
of  F33615-77-C-5221.  The  prime  contract  is  funded  by  the 
United  States  Air  Force  Aero-Propulsion  Laboratory. 

This  report  covers  work  conducted  during  the  period  of 
October  1977  to  November  1979  and  is  part  of  Task  IIIC. 

The  GE  Program  Manager  was  Mr.  Joe  McKenzie  and  the 
Principal  Investigator  was  Mr.  A1  Storace.  The  work  reported 
herein  was  performed  under  the  direction  of  Dr.  Louis  I.  Boehman 
of  the  Mechanical  Engineering  Department  and  Mr.  Antonios  • 

Challita  of  the  Experimental  and  Applied  Mechanics  Division, 
University  of  Dayton  Research  Institute. 

Program  management  for  the  University  was  provided  by 
Mr.  Robert  Bertke. 

This  report  covers  work  conducted  for  project  3066,  task  12, 
entitled  Foreign  Object  Impact  Design  Criteria.  The  contract  was 
sponsered  by  the  Aero  Propulsion  Laboratory,  Air  Force  Systems  Command, 
Wright-Patterson  AFB ,  Ohio  45433  under  the  direction  cf  Sandra  K.  Drake 
( AFWAL/POTA) ,  Project  Engineer. 
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A  MODEL  FOR  PREDICTING  BIRD  AND  ICE 
IMPACT  LOADS  ON  STRUCTURES 


SECTION  1 

INTRODUCTION  AND  SUMMARY 

1.1  BACKGROUND 

Hail,  ice  sheets,  and  birds  are  periodically  ingested  into 
aircraft  engines  during  take-off,  flight,  and  landing  operations. 
The  resulting  damage  to  aircraft  components  such  as  engine  fan 
blades  and  aircraft  windshields  can  lead  to  destruction  of  the 
aircraft  and  crew.  One  of  the  most  serious  threats,  especially 
in  high  speed  flight  at  low  altitudes,  is  bird  ingestion  inf  the 
engine.  The  elements  of  a  jet  engine  which  are  most  vulner?  * 

to  impact  of  ice  and  birds  (foreign  object  damage  (FOD)]  are  ,e 

first-stage  fan  blades.  Conventional  jet  engines  currently 
use  sustain  relatively  minor  damage  with  only  rare  occurrences  of 
catastrophic  failure  when  birds  or  ice  are  ingested.  The  first 
stages  of  these  engines  have  "thick"  titanium  or  stainless  steel 
blades  and  run  at  moderate  rotational  speeds  (chiefly  limited  by 
the  strength-to-weight  ratio  of  the  blade  materials).  However, 
advanced  engines  currently  under  development  require  thin  blades 
with  sharp  leading  edges  and  high  rotational  speeds  in  order  to 
obtain  high  speed  aerodynamic  efficiency.  Efforts  to  increase 
the  performance  of  conventional  engines  envision  the  use  of 
lightweight  composite  materials  to  achieve  higher  rotational  fan 
speeds  and  higher  power-to-weight  ratios.  Both  of  these  direc¬ 
tions  in  fan  jet  engine  evolution  pose  severe  design  problems  for 
the  successful  development  of  a  new  generation  of  FOD-resistant 
fan  blades. 

Progress  in  the  development  of  FOD-resistant  fan  blades 
has  been  hindered  by  the  lack  of  understanding  of  the  mechanisms 
of  FOD  failure,  and  the  lack  of  a  blade  analysis  tool  capable  of 
predicting  the  response  of  a  fan  blade  to  ice  and  bird  impact 
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loading.  The  recent  development  of  sophisticated  dynamic  struc¬ 
tural  analysis  programs  based  on  the  finite-element  methodology 
gives  the  fan  blade  designer  the  basic  tool  needed  to  determine 
structural  response  due  to  impact  loading.  What  is  still  required 
in  order  to  use  and  capitalize  on  the  power  of  these  well- 
developed  finite-element  computer  codes  for  blade  analysis,  is  a 
physically  correct  and  accurate  model  for  the  bird  or  ice  impact 
process. 

This  report  describes  the  development  of  a  bird  and  ice 
impact  loading  model  which  is  specifically  designed  to  be  used 
with  finite-element  structural  analysis  computer  codes.  Two 
recent  developments  from  POD  investigations  have  provided  the 
information  sorely  needed  to  obtain  a  fundamental  understanding 
of  the  FOD  impact  process.  The  first,  and  most  important,  obser¬ 
vation  is  that  bird  and  ice  impacts  are  primarily  fluid  dynamic 
in  nature  and  that  viscous  effects  can  be  ignored.  The  second 
observation,  which  is  particularly  important  for  modelling  impacts 
on  fan  blades,  has  to  do  with  the  cutting  action  of  a  blade  during 
the  slicing  of  the  impacting  object  by  the  leading  edge  of  a  blade 
It  was  observed  that  no  significant  loading  is  attributable  to 
the  cutting  action  itself;  that  is,  the  slicing  force  is  pri¬ 
marily  due  to  the  change  in  direction  which  the  blade  imparts  to 
the  slice  mass.  This  second  observation  basically  means  that  the 
entire  impact  process  of  a  bird,  hail,  or  slab  of  ice  striking  a 
fan  blade  can  be  reasonably  modelled  by  the  methods  of  fluid 
dynamic  analysis  without  having  to  consider  the  tensile  or 
compressive  strength  of  the  impacting  object,^ 

Because  the  model  for  bird  and  ice  impact  loading  described 
in  this  report  is  specifically  designed  to  interface  with  finite 
element  transient  structural  analysis  computer  codes,  the  model 
capitalizes  on  the  fact  that  the  impacted  surface  is  fully 
described  mathematically  in  the  structural  analysis  code  at  any 
instant  of  time.  This  fact,  coupled  with  the  observation  that 
the  impact  process  can  be  modelled  as  an  ideal  ( non-viscous )  fluid 
flowing  onto  the  blade  surface  suggests  that  the  well-developed 


methods  of  ideal  fluid  dynamic  analysis  (potential  flow  theory) 
can  be  used  to  model  FOD  impacts  and  that,  in  particular,  the 
surface  singularity  method  for  solving  complex  potential  flows  is 
the  ideal  tool  to  be  used  for  the  impact  loading  model. 

This  report  describes  how  the  surface  singularity  technique 
is  used  to  compute  the  loads  exerted  on  a  surface  during  FOD 
impacts,  and  how  this  technique  interfaces  with  the  finite-element 
structural  analysis  method.  It  should  be  mentioned  here,  at  the 
outset,  that  the  FOD  loading  model  described  in  this  report  is  not 
based  on  mathematically  exact  solutions  to  a  well  posed  potential 
flow  problem.  Rather,  because  of  computer  time  and  memory  limi¬ 
tations,  only  an  approximate  solution  is  obtained  to  model 
equations.  The  equations  only  approximate  the  true  fluid  dynamic 
event. 

The  overall  model,  while  not  stated  or  solved  with  exact 
mathematical  rigor,  does  include  descriptions  of  the  most  impor¬ 
tant  physical  phenomena  associated  with  FOD  impacts.  These 
physical  phenomena  include  the  following: 

1.  A  true  three-dimensional  treatment  of  the  impact 
process  is  utilized. 

2.  The  shape  and  size  (i.e.,  geometry)  of  the  slice  mass 
is  computed. 

3.  The  impacted  surface  shape  is  arbitrary  and  deformation 
under  load  is  considered  (i.e.,  coupling  between  the 
loading  process  and  target  response  is  treated  in  a 
physically  meaningful  manner). 

4.  Initiation,  duration,  and  termination  of  the  loading 
process  is  described  with  separate  descriptions  used 
for  birds,  ice  spheres,  and  ice  slabs. 

In  the  remainder  of  Section  I  of  this  report,  the  approxima¬ 
tions  used  in  the  formulation  of  the  loading  model  are  discussed 
and  justified.  Section  II  describes  the  methods  and  analysis 
used  to  generate  a  description  of  the  slicing  process  and  how  the 
slicing  process  is  viewed  and  modelled  for  birds,  ice  spheres,  and 
ice  slabs  impacts.  Section  III  of  the  report  describes  the  surface 
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singularity  technique  and  how  it  is  applied  in  the  loading  model* 
Section  IV  is  basically  a  description  of  an  FOD  loading  model 
computer  program,  and  how  it  is  interfaced  with  the  finite-element 
structural  analysis  method. 

1.2  FLUID  DYNAMIC  NATURE  OF  FOD  IMPACTS 

The  impact  of  a  bird  or  a  slab  of  ice  onto  a  jet  engine  fan 
blade  is  a  rather  unusual  impact  problem.  It  involves  a  number 
of  effects  that  are  normally  avoided  in  impact  investigations. 
Firstly,  all  bird  and  ice  impacts  on  jet  engine  fan  blades  are 
oblique.  Obliquity  has  the  effect  of  making  the  problem  three- 
dimensional.  Only  under  very  severe  limitations  ("spherical" 
birds  or  ice  and  nondeforming  targets)  can  the  process  be  reduced 
to  two  dimensions.  The  rigorous  analysis  of  most  truly  three- 
dimensional  impact  events  appears  to  be  beyond  the  current  state 
of  the  art.  The  only  analytic  techniques  which  promise  to  be 
capable  of  addressing  general  oblique  impacts  are  finite  dif¬ 
ference  methods.  At  the  moment  they  are  prohibitively  expensive 
and  only  moderately  accurate.  The  second  important  effect  in 
bird  and  ice  impacts  with  fan  blades  is  that,  in  general,  the 
impact  is  a  slicing,  edge  impact.  This  effect  ensures  three- 
dimensionality,  even  under  the  assumptions  which  render  oblique 
impacts  on  a  plane  surface  two  dimensional.  There  appear  to  be 
no  proven,  existing  analytic  techniques  capable  of  rigorously 
treating  an  oblique,  slicing  edge  impact.  The  problem  must  be 
treated  with  some  degree  of  approximation.  When  compliant 
targets  are  considered,  the  degree  of  approximation  required 
to  obtain  a  solution  increases.  This  section  of  the  report 
describes  an  approach  to  the  solution  of  this  problem  which 
incorporates  rigorously  the  fluid  dynamic  nature  of  the  impact 
event. 

1  2 

It  is  known  from  several  other  studies  '  that  a  material 
like  gelatin  with  10  percent  porosity  adequately  simulates  the 
impact  properties  of  real  birds.  The  early  efforts  in  the  deve¬ 
lopment  of  the  substitute  bird  material  for  use  in  bird-impact 
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studies  were  reported  by  Allcock  and  Collin. 


Wilbeck*  theoretically  and  experimentally  studied  the 
impact  behavior  of  low-strength  projectile  materials  which  he 
characterized  as  soft-body  materials.  These  materials  have  a 
much  lower  strength  than  that  of  typical  target  materials.  An 
impact  involving  a  projectile  of  soft-body  material  against  a 
target  surface  made  of  aluminum  or  steel  generates  stresses  that 
substantially  exceed  the  strength  of  the  projectile  material  but 
are  well  below  the  strength  of  the  target  material. 

Wilbeck*  and  Barber^  have  established  in  their  reports 
that/  unlike  impacts  involving  projectiles  of  strong  materials/ 
the  impact  of  projectiles  of  soft-body  materials  is  dominated  by 
the  tendency  of  the  projectile  material  to  behave  like  a  fluid 
during  the  impact.  Such  soft-body  materials  include/  in  addition 
to  bird  and  substitute  bird  materials/  ice  in  the  form  of  hail  or 
ice  sheets  which  break  off  the  engine  nacelle. 

It  is  useful  to  review  here  the  rationale  for  considering 
bird  and  ice  impacts  as  fluid  dynamic  phenomena.  Recognition  of 
this  type  of  behavior  greatly  aids  in  the  development  of  a 
theoretical  understanding  of  the  impact  process. 


1.2.1  Observations  Based  on  Experimental  Measurements 


Wilbeck  treated  the  impact  of  a  bird  on  a  rigid 
plate  as  an  unsteady  fluid  dynamic  process  and  developed  a 
simplified  one-dimensional  analysis  of  a  homogeneous  right- 
circular  cylinder  of  soft-body  material  impacting  normally  on  a 
rigid  plate.  The  analysis  showed  the  entire  impact  process  to 
occur  in  four  distinct  phases.  In  the  first  phase,  which  is  the 
initial  impact  phase,  very  high  shock,  or  Hugoniot,  pressures  are 
generated.  He  calculated  this  pressure,  using  the  Hugoniot  rela¬ 
tion  for  a  mixture,  together  with  the  shock  properties  of  gelatin. 
The  results  of  the  peak  pressure  measurements  for  the  impact  of 
right  circular  cylinders  of  gelatin  with  10  percent  porosity 
reported  by  Barber,  et  al.^  showed  that  the  measured  pressures 

were  in  good  agreement  with  the  impact  pressures  calculated  by 
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Wilbeck.  More  recent  work  by  Bauer  and  Barber  '  shows  even 
better  agreement.  This  agreement  in  the  initial  phase  of  the 
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impact  process  indicates  that  the  hydrodynamic  description  of  the 
event  is  well  justified. 


After  the  first  phase,  the  high  shock  pressures 
decay  to  steady  fluid  dynamic  pressures.  These  pressures  can  be 
calculated  by  considering  the  process  as  steady  jet  flow.  This 
theoretical  conclusion  is  again  in  agreement  with  measured 
steady-flow  pressures.1'3'5 

Furthermore,  the  calculations  of  shock  decay  by 
Wilbeck1  have  established  that  for  a  projectile  with  a  length- 
to-diameter  ratio  larger  than  a  critical  value,  the  shock  will 
be  severely  weakened  by  radial  expansion  waves  and  the  projec¬ 
tile  should  undergo  complete  shock  decay  to  steady  flow.  Steady 
flow  would  be  expected  to  prevail  if  the  length-to-diameter 
ratio  of  the  projectile  in  the  direction  of  the  impact  exceeds 
approximately  unity.  For  real  birds  striking  end-on,  the 
length-to-diameter  ratio  ranges  from  2  to  3  and  a  steady-flow 
regime  should  occur.  This  conclusion  has  been  amply  supported 
by  the  measurements.1'3'5 

Further  evidence  of  the  tendency  of  the  projectile 
material  to  flow  radially  outward  at  the  impact  location  is 
apparent  during  the  steady-flow  phase.  As  the  radial  pressures 
decrease  during  the  shock  decay,  shear  stresses  develop  in  the 
projectile  material.  If  the  shear  strength  of  the  projectile 
material  is  large  enough  to  withstand  these  shear  stresses,  the 
radial  motion  of  the  projectile  will  be  impeded.  On  the  other 
hand,  the  projectile  material  will  begin  to  flow  if  its  shear 
strength  is  smaller  than  the  shear  stresses  developed.  The 
experiments  have  confirmed  that  for  real  birds,  gelatin  and  ice, 
the  shear  strength  is  low  enough  for  the  pressures  generated  to 
cause  the  projectile  material  to  flow. 

1.2.2  Simplified  Potential  Flow  Model  and  Comparisons 

Of  Predicted  and  Measured  steady-state  pressure 

Distributions 

The  absence  of  any  entrainment  of  the  surrounding 
fluid  (air)  during  the  impact  process  led  to  the  recognition 


that  the  steady-flow  phase  o£  the  bird  impact  would  be  ideally 
suited  for  modeling  by  potential  flow  theory.  This,  in  turn,  led 
to  the  development  of  a  three-dimensional  potential  flow  model 3 
for  predicting  the  pressure  distribution  produced  by  the  steady 
flow  of  a  cylindrical  jet  impacting  obliquely  on  a  flat  plate. 

The  solution  of  the  three-dimensional  Laplace's  equation  for  a 
steady,  incompressible,  irrotational  flow  was  obtained  by  the 
superposition  of  the  two  elementary  solutions.  These  elemen¬ 
tary  solutions  were  the  uniform  flow  of  a  fluid  in  a  round  duct, 
and  the  uniform  distribution  of  planar  sources  over  the  elliptical 
area  on  the  target  defined  by  the  intersections  of  the  round  duct 
and  the  flat  plate.  The  details  of  this  procedure,  which  is  an 
elementary  application  of  the  surface  singularity  technique,  are 
described  in  Reference  3 . 

The  predicted  steady-flow  pressure  distribution 
calculated  showed  very  good  agreement  with  the  measured  values 
over  the  central  portions  of  the  jet  flow.3  However,  the 
agreement  was  rather  poor  near  the  edges  of  the  jet,  and  this 
descrepancy  was  essentially  due  to  the  fact  that  the  vorticity 
effects  which  are  important  at  the  edges  of  the  jet  are  not 
properly  modelled  by  the  simple  superposition  of  two  elementary 
potential  flow  solutions.  Nevertheless,  the  fact  that  this 
simple  model  correctly  predicted  pressure  distributions  over 
the  major,  and  most  important,  region  of  the  impact  zone  for  a 
wide  range  of  impact  angles  further  supports  the  concept  that 
bird  and  ice  impacts  are  essentially  ideal  fluid  flow  in  nature. 

1.3  ASSUMPTIONS  AND  THEIR  JUSTIFICATIONS 

It  is  clear  from  the  foregoing  that  the  impact  of  soft-body 
materials  such  as  birds,  realistic  substitute  bird  materials  like 
gelatin,  and  ice  at  the  typical  velocities  of  impact  may  be 
analyzed  within  the  framework  of  hydrodynamic  theory.  The 
impact  process  is  modelled  in  terms  of  normal  and  oblique 
impact  of  a  right  circular  cylinder  against  rigid  and  compliant 
target  surfaces.  The  theoretical  analysis,  of  course,  requires 
several  simplifying  assumptions. 
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The  projectile  material  is  considered  to  be  homogeneous 
even  when  large  amounts  o£  porosity  are  considered.  This 
assumption  would  appear  to  be  somewhat  unrealistic  in  the  case 
o f  real  birds,  although  extensive  test  data  do  not  indicate 
gross  inhomogeneties The  assumption  is  quite  reasonable  for 
the  substitute  material  (gelatin  with  10  percent  microballoon)  and 
ice. 

In  the  analysis  of  the  impact  process ,  the  strength  of 
the  projectile  material  is  considered  to  be  negligible.  This 
assumption  is  reasonable  for  the  typical  projectile  and  target 
materials  of  interest  at  the  impact  velocities  of  interest. 

It  is  assumed  that,  at  least  as  a  first  approximation,  the 
fluid  flow  may  be  treated  as  incompressible.  This  is  a  reasonable 
assumption  in  view  of  the  fact  that  the  measured  steady-state 
pressures  are  quite  small  in  comparison  to  the  pressures  required 
to  produce  significant  density  changes  in  water. 

In  modeling  the  steady-flow  phase  by  potential  flow  anal¬ 
ysis,  it  is  assumed  that  the  flow  is  inviscid  and  irrotational . 

Since  the  impact  process  does  not  entail  any  entrainment  of  the 
surrounding  fluid,  and  since  the  time  over  which  steady  flow 
exists  is  small  compared  to  the  time  required  to  establish  strong 
vorticity  in  the  flow,  it  seems  reasonable  to  treat  the  flow 
field  as  effectively  irrotational,  at  least  over  the  central 
portion  of  the  jet.  However,  the  fact  that  this  simplification 
fails  to  model  adequately  the  edges  of  the  jet  is  well-recognized. 

It  must  be  noted  that  the  presence  of  a  significant  amount 
of  porosity  in  birds  results  in  a  very  low  sonic  velocity  (of  the 
order  of  40  m/s  for  gelatin  with  10  percent  porosity).  Then,  for 
the  typical  velocities  at  which  bird  impact  occurs,  the  initial 
impact  process  will  be  largely  supersonic.  Thus,  while  the 
shock  wave  is  weakened  by  the  release  waves  and  is  ultimately 
eliminated  for  a  subsonic  impact,  it  will  not  disappear  for  a 
supersonic  impact.  The  shock  propagation  velocity  will  decrease 
until  it  becomes  equal  to  the  impact  velocity,  then  a  standing  shock 
will  be  established.  Behind  this  shock  the  flow  will  be  subsonic 
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and  will  follow  steady  flow  streamlines.  Wilbeck's  calculations1, 
based  upon  potential  flow  theory  for  a  supersonic  bird  impact  at 
normal  incidence,  show  that  the  steady-flow  pressure  at  the 
center  of  impact  is  almost  independent  of  porosity,  suggesting 
that  the  decrease  in  density  due  to  porosity  is  offset  by  the 
increase  in  compressibility. 

1.4  SUMMARY  DESCRIPTION  OF  LOADING  MODEL 

If  the  target  deforms  locally  during  the  impact,  as  well 
as  translating  and  rotating,  then  a  loading  model  capable  of 
generating  the  local  distribution  of  loading  pressures  during 
the  impact  is  required.  As  described  previously,  birds  behave 
like  a  fluid  during  impact  and  the  distribution  of  surface 
pressure  during  impact  is  directly  related  to  the  fluid  nature 
of  the  event.  To  successfully  predict  the  surface  pressure 
distributions  on  deforming  surfaces,  it  is  necessary  to  utilize 
a  fundamental  fluid  dynamic  approach.  The  shear  stress  distribu¬ 
tions,  due  to  boundary  layer  {viscous)  effects  can  safely  be 
ignored  for  both  normal  and  oblique  impacts  since  the  ratio  of 
maximum  normal  stress  (pressure)  to  maximum  shear  stress  is  of 
the  order  of  one  over  the  skin  friction  coefficient.  Therefore, 
the  problem  is  reduced  to  determining  the  pressure  distribution 
over  the  impact  area. 

If  the  event  is  assumed  to  be  dominated  by  the  quasi-steady 
flow  portions  of  the  impact,  as  described  previously,  the  process 
can  be  thought  of  as  the  flow  of  a  jet  onto  a  surface  as  illustrated 
in  Figure  1.  (Shock  effects,  if  demonstrated  to  be  important,  can 
be  evaluated  separately  and  superimposed  on  the  quasi-steady 
flow  results.) 

The  characteristic  (and  maximum)  pressure  in  quasi-steady 
fluid  flow  is  the  Bernoullian  stagnation  pressure  (1/2  pvr2) 
and  the  important  independent  parameters  are  the  impactor  density, 

0,  and  the  impact  velocity,  Vr.  The  loads  are  specified  if  the 
stagnation  pressure  and  pressure  coefficient  distribution  can  be 
determined.  As  pointed  out  previously,  porosity  (of  birds)  has 
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MAJOR  AXIS 


Figure  1.  Oblique  Impact 

little  effect  on  the  stagnation  pressure.  The  zero  porosity 
density  (Of)  and  incompressible  flow  assumptions  can,  therefore, 
be  applied.  Thus,  the  problem  is  reduced  to  the  classical 
hydrodynamic  problem  of  the  flow  of  a  liquid  jet  onto  a  surface. 
The  problem  can  be  further  simplified  by  making  the  assumption 
of  irrotational  flow:  an  assumption  which  is  supported  by  the 
fact  that  in  the  impact  region  fluid  inertia  forces  dominate  over 
viscous  forces.  With  the  assumptions  of  steady,  incompressible, 
irrotational  flow,  the  problem  may  be  treated  as  a  steady,  poten¬ 
tial  flow  described  by  Laplace's  equation.  A  serious  complication 
remains  however.  The  boundary  of  the  jet  is  a  free  streamline 
whose  position  is  not  known  "a  priori.”  The  static  pressure  is 
continuous  across  a  free  streamline  while  the  velocity  tangent 
to  this  streamline  and  the  stagnation  pressure  are  both  discon¬ 
tinuous.  Along  a  free  streamline,  the  magnitude  of  the  velocity 
is  constant  according  to  Bernoulli's  equation.  Thus,  even  though 


the  governing  partial  differential  equation  (Laplace's  equation) 
is  linear,  the  free  streamline  boundary  condition  is  nonlinear. 
These  conditions  can  be  described,  in  two-dimensional  flows,  with 
complex  variable  theory  and  exact  analytical  solutions,  and  have 
been  obtained,  for  two-dimensional  jets  impacting  on  flat  surfaces. 
However,  no  three-dimensional  jet  impact  flows  (such  as  those  of 
interest  here)  have  been  solved  analytically.  Not  even  the  case 
of  an  axisymmetric  (circular)  jet  impacting  normal  to  a  flat 
plate  has  been  solved,  because  complex  variable  theory  becomes 
extremely  cumbersome  in  cylindrical  coordinate  systems.^  Approx¬ 
imate  analytical  solutions  are  available  for  predicting  the 

pressure  distribution  for  circular  jets  impacting  normal  to  flat 
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plates.  '  Reference  8  contains  semi-empirical  expressions  and 
techniques  for  computing  pressure  distributions  for  cylindrical 
jets  impacting  obliquely  on  flat  surfaces.  The  general  problem 
of  a  jet  of  an  arbitrary  cross  section  impacting  on  an  arbitrary 
surface  can  only  be  solved  by  numerical  methods. 

Numerical  potential  flow  solutions  for  jets  (jets  bounded 
by  free  streamlines)  are  very  difficult  to  obtain  because  the 
location  of  the  free  streamlines  are  not,  in  general,  known. 

For  circular  jets,  at  normal  incidence,  the  position  of  the  free 
streamlines  can  be  adequately  assumed  in  order  to  start  a 
solution.  For  oblique  impacts,  particularly  on  curved  surfaces, 
assuming  the  free  streamline  position  would  be  a  major  under¬ 
taking  in  addition  to  the  considerable  task  of  numerically 
solving  Laplace's  equation  in  three  dimensions.  Accordingly, 
an  approximate  numerical  solution  which  would  provide  a  reason¬ 
ably  accurate  description  of  the  pressure  distribution  (but  not 
necessarily  satisfying  the  free  streamline  condition)  was 
developed . 

Three-dimensional  potential  flow  theory  was  used  to  develop 
a  model  for  predicting  the  pressure  distribution  produced  by  the 
steady  flow  of  a  cylindrical  jet  impacting  on  a  flat  plate.  It 
was  assumed  that  pressure  distribution,  as  calculated  for  this 
fluid  dynamic  problem,  would  provide  a  reasonable  description 
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of  the  steady  flow  portion  of  a  bird  impact.  Incompressible, 
irrotational  flow  assumptions  were  made,  but  the  free  streamline 
boundary  condition  was  not  imposed.  The  method  of  superposition 
of  solutions  was  used  even  though  superposition  is  not  strictly 
valid  (although  Laplace's  equation  is  linear,  the  free  stream 
boundary  condition  is  nonlinear) .  The  numerical  approach  used 
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was  based  on  the  method  of  surface  singularities. 

The  approximate  numerical  solution  involved  the  super¬ 
position  of  two  elementary  solutions  to  Laplace's  equation. 

The  two  solutions  used  were:  (1)  the  uniform  flow  of  a  fluid 
in  an  infinite  jet  of  arbitrary  but  constant  cross-sectional 
area;  and  (2)  the  uniform  distribution  of  planar  fluid  sources 
over  the  impact  area.  The  velocity  of  the  fluid  within  the 
infinite  jet  was  considered  constant,  while  the  fluid  outside 
this  region  was  assumed  to  be  at  rest.  Surface  sources  were 
distributed  on  the  target  surface.  The  strength  of  the  surface 
sources  was  selected  to  provide  the  correct  boundary  condition 
at  the  target  (no  flow  normal  to  the  target  surface).  The 
complete  solution  to  this  problem  for  jets  of  circular  cross 
sections  impacting  obliquely  on  rigid  flat  surfaces  was  devel¬ 
oped  by  Boehman  of  the  University  of  Dayton  and  is  presented 
in  Reference  3. 

This  model  can  be  used  for  any  arbitrarily  shaped  impact 
surface.  The  impact  area  is  divided  into  small  flat  elements, 
and  a  uniform  distribution  of  sources  is  assumed  to  cover  each 
element.  During  the  analysis  of  an  impact,  in  which  local  target 
deformation  takes  place,  the  deformed  target  shape  in  the  impact 
zone  is  calculated  at  each  time  interval  employed  in  the  dynamic 
structural  analysis.  The  geometry  of  the  impact  zone  can  then 
be  provided  at  each  time  interval  to  the  loading  model,  and  the 
pressure  distribution  appropriate  to  the  target  at  that  time  can 
be  determined.  As  the  structural  analysis  calculation  proceeds, 
the  local  shape,  location,  and  velocity  of  the  impact  area  is  up¬ 
dated  and  made  available  to  the  loading  model.  The  loading  model, 
in  turn,  provides  updated  pressure  and  pressure  distribution 
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information  for  the  structural  response  computation.  The  loading 
model  is,  thus,  fully  interactive  with  the  structural  response 
calculation.  The  duration  of  the  impact  is  computed  by  computing 
the  amount  of  slice  consumed  as  a  function  of  time  and  setting 
surface  pressures  equal  to  zero  after  the  slice  is  consumed. 

The  loading  model  is  capable  of  detailed  interaction  with 
the  structural  response  model  and  is  capable  of  dealing  with 
target  translation,  rotation,  and  local  deformation.  The  load/ 
response  coupling  modelled  in  this  formulation  should  be  capable 
of  accurate  prediction  of  pressure  for  both  overall  target 
response  and  local  deformation. 

The  principal  limitation  of  the  loading  model  is  that 
it  does  not  include  transient  effects.  The  most  significant 
transient  effect  is  the  shock  process.  However,  the  porosity 
present  in  birds  appreciably  reduces  the  shock  stresses,  while 
the  steady  flow  pressures  are  not  significantly  affected.  In 
addition,  impact  obliquity  (bird-blade  impacts  are,  in  general, 
oblique)  also  reduces  the  relative  importance  of  shock  stresses. 
Therefore,  it  is  not  obvious  that  neglect  of  the  shock  aspects 
of  bird  impact  on  blades  is  a  significant  deficiency.  Another 
transient  aspect  of  ice  sphere  impacts  is  the  variation  of  a 
slice  cross-sectional  area  at  the  target  surface  during  impact. 
This  results  mainly  in  a  time  variation  of  the  impact  area.  If 
the  flow  remains  quasi-steady  during  these  variations  (i.e., 
the  "velocity*  of  the  variation  is  low  with  respect  to  local 
sound  speed;  probably  a  good  assumption),  then  the  model  can  be 
modified  to  describe  these  effects.  The  size  and  geometry  of 
the  "jet"  which  flows  onto  the  target  surface  must  be  updated 
in  time  incrementally  to  describe  the  variation  of  impact  area 
with  time. 

1 . 5  SUMMARY 

In  the  current  version  of  the  loading  model,  the  pressure 
distribution  is  based  on  a  steady  flow  analysis.  This  method 
of  analysis  is  extended  into  a  dynamic  analysis  for  deforming 
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targets  by  treating  the  flow  at  any  instant  of  time  as  a  quasi¬ 
steady  flow  using  the  instantaneous  relative  velocity  between 
the  bird  slice  and  the  deforming  target  as  the  characteristic 
velocity.  The  instantaneous  shape  of  the  target  is  used  to 
define  the  surface  on  which  impact  occurs. 

In  the  analysis  of  the  impact  process,  the  strength  of 
the  projectile  material  is  considered  to  be  negligible.  This 
assumption  is  quite  reasonable  for  the  typical  projectile  and 
target  materials  of  interest. 

It  is  assumed  that,  at  least  as  a  first  approximation,  the 
fluid  flow  may  be  treated  to  be  incompressible.  This  is  a 
reasonable  assumption  in  view  of  the  fact  that  the  measured 
steady-state  pressures  are  quite  small  in  comparison  to  the 
pressures  required  to  produce  significant  density  changes  in 
birds  or  ice.* 

By  virtue  of  the  incompressible  flow  and  steady  flow 
assumptions,  together  with  negligible  strength  of  the  projectile 
material,  the  problem  of  predicting  soft  body  impact  loads  is 
amenable  to  analysis  within  the  following  real  constraints. 

The  time  required  to  obtain  a  reasonable  solution  must  be  small 
in  comparison  to  the  time  required  to  compute  the  structural 
response,  and  the  computer  storage  requirement  for  the  loa  ing 
model  must  be  small  in  comparison  to  the  structural  analysis 
computer  program  storage  requirements. 

Three-dimensional  potential  flow  theory  was  chosen  as  the 
most  appropriate  method  for  modelling  the  impact  process.  Our 
initial  work  in  using  this  approach  was  quite  successful. ^  In 
Reference  3,  the  method  of  surface  singularities  was  used  to 
determine  the  velocity  and  pressure  fields  due  to  circular  jets 
impacting  at  oblique  angles  on  flat  rigid  plates.  This  same 
technique  was  used  to  develop  the  loading  model  described  in 
this  report.  The  method  developed  in  Reference  3  was  generalized 
to  include  impact  surface  curvature  and  deformation  velocity, 
arbitrary  cross-sectional  area  of  the  impacting  flow,  and  a 
method  for  generating  planar  quadrilateral  surface  elements 
given  finite-element  surface  nodal  point  locations. 
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SECTION  2 


DISCUSSION  OF  BIRD  AND  ICE  IMPACT 
ON  MOVING  AND  STATIONARY  BLADES 

While  the  discussion  in  the  previous  section  dealt  largely 
with  the  modeling  o£  the  soft-body  impacts  on  surfaces  (bird- 
strike  on  aircraft  windshield/  for  example)/  the  other  area  of 
major  concern  in  the  study  of  these  impacts  is  the  ingestion  of 
foreign  bodies  in  jet  engines  resulting  in  impacts  on  fan  or 
compressor  blading.  Typically,  this  problem  might  arise  from 
the  ingestion  of  a  bird  or  of  ice  sheets  breaking  off  the 
nacelle.  The  primary  consequences  of  such  impact  would  be 
deflection,  bending,  and  rupture  of  the  impacted  blades.  Although 
this  can  give  rise  to  secondary  impacts  and  related  effects,  the 
present  research  program  is  confined  only  to  the  primary  foreign 
object/blade  impacts.  It  is  useful  to  point  out  that  while  the 
subsequent  discussion  might  explicitly  deal  with  the  bird/blade 
impact  process,  the  analysis  and  the  conclusions  herein  are  not 
restricted  to  bird  impact  alone  but  apply  to  any  foreign  body 
impact  (it  is  assumed,  of  course,  that  we  would  still  be  con¬ 
sidering  only  "soft-body"  materials) . 

The  bird/blade  impact  process  differs  in  an  essential  way 
from  the  earlier-discussed  bird  impact  over  extended  surfaces. 
Unlike  the  bird  striking  an  aircraft  windshield,  for  instance,  a 
bird  entering  an  engine  is  cut  into  slices  by  the  chopping  action 
of  the  first-stage  fan  blades.  Thus,  before  one  proceeds  with 
the  analysis  of  the  impact  process,  it  is  necessary  to  first 
establish  the  geometry  and  parameters  of  the  slices  of  bird,  or 
ice,  formed  by  a  rotating  fan  stage. 

2.1  SLICING  MODEL  DEVELOPMENT 

During  1976  a  few  experiments  were  conducted  at  the  Air 
Force  Materials  Laboratory  by  the  University  of  Dayton  Research 
Institute  to  investigate  the  effects  of  slicing  edge  impacts 
on  the  impact  loads.  Birds  were  fired  over  thin  wedges  and 
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cylindrical  wires  mounted  in  a  ballistic  pendulum.  From 
measurements  of  the  momentum  transfer  during  the  slicing  process, 
estimates  of  the  slicing  force  could  be  made.  It  was  found  that 
the  slicing  force  was  predominantly  due  to  fluid  dynamic  drag 
with  no  significant  loading  attributable  to  the  cutting  action 
itself.  These  experiments  suggest  that  it  is  reasonable  to 
ignore  the  cutting  forces  in  comparison  to  the  forces  required 
to  decelerate  and/or  change  the  direction  of  an  impacting  soft 
body.  Therefore,  the  slicing  model  development  reduces  to 
simply  a  geometric  problem  of  determining  the  dimensions  and 
weight  of  a  slice.  For  birds  and  ice  spheres  (such  as  hail) 
the  velocity  of  the  ingested  object  is  small  compared  to  the 
velocity  of  the  aircraft;  thus,  the  object  velocity  can  be 
ignored.  For  ice  slabs,  such  as  ice  breaking  loose  from  an 
engine  nacelle,  the  velocity  of  the  slab  relative  to  the  nacelle 
is  not  well  defined.  A  reasonable  assumption,  at  least  for  short 
inlet  systems,  is  to  ignore  the  velocity  of  the  slab  relative  to 
the  nacelle.  For  long  inlet  systems,  it  might  be  possible  to 
estimate  the  relative  velocity  between  slab  and  nacelles  by 
considering  the  action  of  aerodynamic  drag  forces  on  an  ice 
slab  after  it  has  broken  loose  from  a  forward  section  of  a 
nacelle . 

2.1.1  Slicing  Model  Development  for  Birds 

A  bird  is  idealized  as  a  right,  circular  cylinder 
with  a  length-to-diameter  ratio  of  2  and  the  velocity  of  the 
bird,  relative  to  the  aircraft,  is  taken  to  be  equal  and  opposite 
to  the  aircraft  velocity. 

In  the  following  analysis,  a  coordinate  system 
attached  to  the  blade  is  used.  The  bird/blade  interaction 
geometry,  in  such  a  coordinate  system,  is  shown  in  Figure  2.  The 
following  information  is  assumed  to  be  known  (supplied  as  input 
data  to  the  loading  model  computer  program); 

M  -  number  of  blades  per  stage 

n  -  blade  rotational  speed  (rpm) 
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Z.  -  distance  £rom  the  rotational  axis  of  the 
rotor  to  the  point  on  the  blade  at  which 
the  center  of  the  impact  occurs 

<5  -  blade  orientation  angle  (defined  as  the 

angle  5  in  Figure  1) 

Vb  -  axial  velocity  of  the  bird 
-  bird  weight 

Pb  -  bird  density. 

An  infinite  number  of  slice  shapes  is  possible 
for  a  given  set  of  these  input  parameters;  the  shape  depends  on 
the  orientation  of  the  bird  relative  to  the  blade,  and  on  the 
span  location  on  the  blade  at  which  the  impact  occurs.  Since 
worst  case  slice  shapes  are  desired  (i.e.,  slice  shapes  having 
the  largest  possible  slice  mass),  the  orientation  of  the  axis 
of  the  bird  and  the  center  of  impact  of  a  slice  are  chosen  to 
produce  a  slice  having  the  largest  possible  mass.  The  largest 
slice  mass  occurs  when  the  axis  of  the  slice  is  coincident  with 
the  axis  of  the  right  circular  cylinder.  The  slice  geometry 
depicted  in  Figure  2  corresponds  to  this  worst  case  situation. 

With  the  bird  idealized  as  a  right  circular 
cylinder  with  L  ■  2D,  the  diameter  of  the  bird  is  determined 
from  the  bird  weight  and  density  as 


The  tangential  velocity  of  the  bird  Vfc  is  computed 
at  the  impact  radius  Z^  and  is  given  by 


2TrnZ . 


vt  *  60 

and  the  magnitude  of  the  velocity  of  the  bird  relative  to  the 
blade  Vr  is  given  by 

Vr  -j/ 


(2) 


Avb2  +  vt2) 


(3) 
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The  direction  of  Vr  relative  to  the  plane  of 
rotation  of  the  blade  8  at  the  impact  location  Z ^  is 

6  -  sin“1(Vb/Vr)  .  (4) 

Vr  is  assumed  to  have  no  component  in  the  radial  direction.  The 
impact  angle,  defined  as  8  in  Figure  2  is 

0  -  6  -  S  (5) 

and  the  blade  spacing  (defined  as  Z^  and  denoted  by  S  in 
Figure  2)  is  given  by 

S  ■  2irZi/N  .  (6) 

The  bird  slice  width,  h,  is  then  given  by 

2irZ.  V. 

h  ■  Ssin  B  ■  -t? —  rr-  (7) 

r 


Referring  to  Figure  2  we  find  that  the  bird-slice  weight  is 
given  by 

+h/2 

-  v  /  'W77 

-h/2 


The  use  of  Equations  2  and  3  in  Equation  7  shows 
that  the  bird-slice  width  at  any  impact  radius  depends  only  on 
the  blade  parameters,  N  and  n,  and  on  the  aircraft  speed.  The 
bird-slice  weight  at  any  impact  radius  depends,  in  addition  to 
the  above  three  parameters,  only  on  the  density  and  total 
weight  of  the  bird. 


If  h/sin  8  is  greater  than  the  cord  of  the  blade 
at  Z^,  then  not  all  of  the  slice  will  impact  on  the  blade. 

2.1.2  Slicing  Model  for  Ice  Spheres 

For  ice  spheres  (such  as  hail)  the  axial  velocity 
of  the  sphere  relative  to  the  blade  is  taken  to  be  equal  and 
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opposite  to  the  aircraft  velocity,  just  as  for  birds.  The  same 
expressions  for  slice  width  impingement  angle,  and  V„  applies 

9 

as  for  birds.  The  slice  weight  W  is  given  by 


W« 


'is 


sis 


h/2 


2irp 


i  f[ii)2  -  *2]d« 


or 


where 


w 


8is 


(9) 


h  .  D 

7-1 


and  is  the  ice  density.  For  ice  spheres  it  is  presumed  that 
the  diameter  of  the  sphere  is  a  specified  quantity  along  with 
the  ice  density. 

2.1.3  Slicing  Model  for  Ice  Slab 

For  slabs  of  ice  the  operator  specified  information 
is  presumed  to  be  the  length  of  the  slab,  L,  the  thickness  of  the 
slab,  &Z,  (i.e.,  how  much  of  the  blade  span  is  to  be  exposed  to 
the  slab)  and  the  axial  velocity  of  the  slab  relative  to  the 
nacelle  (this  replaces  Vb  in  Equations  3/  4,  and  7).  The  slice 
width,  h,  is  again  given  by  Equation  7  with  replaced  by  the 
axial  velocity  of  the  slab  relative  to  the  nacelle.  The  slice 
mass  is  then  given  by 


W 


Sg  -  Pith)  (A  Z )  L 


(10) 


2.2  BIRD  SLICE  GEOMETRY  PARAMETRIC  STUDY  RESULTS 

Using  the  equations  developed  in  paragraph  2.1.2,  the  bird- 
slice  parameters  were  computed  for  both  starling  impacts  (3-ounce 
birds)  and  big  bird  impacts  (1.5-pound  birds).  Three  different 
blade  configurations,  viz.,  J-79,  APSI,  and  F-101,  are  considered 
and  typical  values  of  the  rotor  speed,  number  of  blades  per  stage, 
blade  orientation  angle,  and  impact  radius  are  used.  The  input 
parameters  and  the  computed  quantities  are  shown  for  the  two  cases 


of  bird  weights  in  Tables  I  and  II.  The  primary  observation  from 
this  parametric  study  is  that  bird  (and  ice  sphere)  impacts  are 
highly  oblique. 


2.3  DURATION  OF  IMPACT 

The  amount  of  a  slice  consumed  at  any  time  during  impact 
is  directly  related  to  the  velocity  of  the  slice  relative  to 
the  blade  (Vr  in  Figure  1). 
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TABLE  I.  BLADE  CALCULATION  INPUT  DATA  SHEET  ( 3-02  BIRD) 
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TABLE  II.  BLADE  CALCULATION  INPUT  DATA  SHEET  (1.5-LB  BIRD) 
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SECTION  3 

POTENTIAL  FLOW  MODELING  OF  FOD  IMPACTS 

The  impact  o f  a  soft  body  material  (bird  or  ice  slices)  on 

a  single  blade  is  modelled  using  a  well  developed  method  of 

potential  flow  theory  known  as  the  surface  singularity  method. 

The  procedure  used  to  apply  the  method  of  surface  singularities 
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to  FOD  impact  is  based  largely  on  the  work  of  Hess  and  Smith. 

3.1  GENERAL  DESCRIPTION  OF  THE  SURFACE  SINGULARITY  METHOD 

The  method  of  surface  singularities  has  its  origin  in  a 
well-known  theorem  of  potential  flow  theory  which ,  in  essence, 
states  that  if  the  velocity  potential  or  its  derivative  is 
known  over  the  entire  boundary  of  a  potential  flow,  then  the 
velocity  field  may  be  determined  throughout  the  region  of  the 
flow. 6  With  this  theorem,  and  with  the  aid  of  Green's  theorem. 
Green's  equivalent  stratum  of  sources  and  doublets  (singularities) 
is  developed.  This  stratum  essentially  states  that  the  velocity 
potential  of  a  fluid  in  motion  can  be  expressed  in  terms  of 
either  a  distribution  of  sources  and/or  doublets  over  the  flow 
boundary.6  For  nonlifting  potential  flows,  only  source  distri¬ 
butions  need  be  considered  whereas  doublet  distributions  are 
required  for  analysis  of  lifting  potential  flows. Source 
distributions  are  considerably  easier  to  implement  and,  based 
upon  the  good  results  obtained  in  our  preliminary  work  using 
source  distributions,^  only  source  distributions  are  used  in 
the  model  described  in  this  report.  The  method  of  surface 
singularities  as  used  in  the  FOD  impact  model  utilizes  a  distri¬ 
bution  of  source  density  on  the  boundaries  of  the  flow,  and 
solves  for  the  distribution  of  source  density  necessary  to 
satisfy  specified  boundary  conditions.  Once  the  source  density 
distribution  is  known,  the  flow  velocities  on  the  boundaries  of 
the  flow  field  and  throughout  the  flow  field  can  be  computed. 

In  order  to  understand  the  method  of  surface  singularities, 
consider  a  steady  flow  of  a  perfect  fluid  impinging  on  a  three- 
dimensional  body  as  shown  in  Figure  3.  in  the  initial  discussion 
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Figure  3.  A  Three-Dimensional  Flow  Field 

o £  the  method,  the  steady  flow  will  be  considered  to  be  infinite 
in  extent.  The  modifications  required  to  adapt  the  method  to 
flows  of  finite  extent  such  as  jet  flows  will  be  discussed  later. 

3.1.1  Surface  Singularity  for  Flows  of  Infinite  Extent 

Let  S  be  the  surface  of  the  body  and  have  an 
equation  of  the  form 


F(x,  y,  *)  ■  0 


(11) 


where  x,  y ,  and  z  are  the  Cartesian  coordinates.  For  simplicity 
assume  that  the  onset  flow  (defined  as  the  flow  field)  is  a  uni¬ 
form  stream  of  unit  magnitude,  and  let  it  be  denoted  by  the  constant 
vector  ^.with  components  V,^,  v-y,  and  V.z,  respectively,  along 
the  coordinates  axes  x,  y,  and  z  where; 
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The-  fluid  velocity  at  a  point  may  ba  axpraasad  as 
the  negative  gradient  of  a  potential  function  ♦  .  Which  satisfies 
Laplace's  Equation  in  the  region  R*  exterior  to  the  surface  S, 
has  a  zero  normal  derivative  on  S,  and  approaches  the  proper 
uniform  stream  potential  at  infinity.  Making  use  of  the  super¬ 
position  principle,  the  potential  function  ♦  is  viewed  as 

*  -  4.  +  4  (13) 

where  is  the  uniform  stream  potential  and  is  equal  to 

4.  "  -<V.xx  +  V„yy  +  V-Xz)  (14) 

and  $  is  the  disturbance  potential  due  to  the  body.  Then  <p 
should  satisfy 

A*  -  0  in  R'  (15) 

■  it  •  grad#  j  ■  n  •  (16) 

■  '  P -0  ■  P-0 

0  for  x2  +  +  z2+  pm  (17) 

where  A  denotes  the  Laplacian  operator  and  ft  is  the  unit  outward 
normal  vector  at  a  point  of  the  surface. 

Prom  potential  theory/  it  can  be  shown  that  this 
potential  may  be  evaluated  in  terms  of  a  surface  source  density 
distribution  with  which  the  body  surface  may  be  considered  to 
be  covered.6  Then  4  may  be  written  as 

+  <x,y,z)  dM  (18) 

where  r(p/q)  is  the  distance  from  the  integration  point/  q,  on 
the  surface  to  the  field  point,  p,  where  the  potential  is  being 
evaluated  as  shown  in  Figure  4.  The  function  o  must  be  deter¬ 
mined  so  that  <p  satisfies  the  normal  derivative  condition. 
Equation  16.  On  the  body  surface  the  normal  derivative  of  $  is 

-  -2s-<r(P) 

s 


^,^7rhy<r(q,‘“  ■ 
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Figure  4.  Notation  used  in  Describing  the  Potential  Due  to  a 
Surface  Source  Density  Distribution 

Substituting  the  value  of  •!£-)  in  Equation  16  gives  the  integral 

o  A  S 

equation  for  a  as 

2TO-(p)  -J),  A(r(p^qy)°’lq)d3  *  ~n(pl'^  <2°> 

This  is  seen  to  be  a  two-dimensional  Fredholm  integral  equation 
of  the  second  kind  over  the  surface ,  S.  Once  this  equation  is 
solved  for  a ,  the  disturbance  potential  $  may  be  evaluated  from 
Equation  18  and  the  disturbance  flow  velocities  from  the  deriva¬ 
tives  of  Equation  19  in  the  coordinate  directions. 

Some  advantages  of  this  method  are:  the  equation 
that  must  be  solved  is  a  two-dimensional  one  over  the  body  sur¬ 
face  rather  than  a  three-dimensional  one  over  the  entire  exterior 
flow  field,  and  the  method  can  be  used  to  calculate  flows  about 
arbitrary  bodies.  There  is  no  restriction  that  the  body  be 
slender,  analytic,  or  simply,  connected  or  that  the  disturbance 
velocities  due  to  the  body  be  small  compared  with  the  velocity 
of  the  onset  flow.  However,  it  is  required  that  the  body  sur¬ 
face  have  a  continuous  normal  vector  n. 


3.1.2  Complications  of  Flow  of  Finite  Extent 


For  flows  of  finite  extent ,  such  as  jet  flows, 
obtaining  an  exact  solution  requires  that  the  boundary  stream¬ 
lines  of  the  flow  must  be  included  as  part  of  the  flow  surface, 

S,  over  which  sources  are  distributed.  The  addition  of  this 
surface  poses  severe  complications.  As  was  mentioned  in  para¬ 
graph  1.4 ,  the  position  of  the  boundary  streamline  of  the  jet  is 
unknown.  Furthermore,  the  jet  surface  streamlines  are  a  special 
type  of  streamline  known  as  free  streamlines  along  which  the 
magnitudes  of  the  velocity  is  constant.  The  difficulties  which 
the  presence  of  a  free  streamline  causes  in  attempting  to  obtain 
exact  numerical  or  analytical  solutions  were  discussed  in  para¬ 
graph  1.4.  Because  of  the  very  real  computer  time  and  computer 
storage  limitations  that  had  to  be  considered  in  the  development 
of  the  present  model,  it  was  not  feasible  to  attempt  an  exact 
numerical  solution.  Fortunately,  reasonably  accurate  pressure 
distributions  can  be  obtained  with  the  surface  singularity  method 
without  attempting  to  satisfy  the  free  streamline  boundary  con¬ 
dition  exactly.  This  is  accomplished  by  use  of  a  suitable  model 
for  the  onset  flow.  Thus,  in  the  F00  loading  model  discussed  in 
this  report,  sources  are  distributed  only  over  the  surface  on 
which  the  impact  occurs.  The  model  thus  utilizes  the  same  basic 
procedure  as  that  described  in  paragraph  3.1.1  for  flows'  of 
infinite  extent,  except  that  the  onset  flow  velocity  distribution 
is  now  not  a  constant  velocity  field  but  is  a  spatially  varying 
velocity  field.  As  is  shown  in  Reference  3,  even  the  simplest 
possible  onset  flow  model,  which  has  a  uniform  distribution  of 
velocity  over  the  cross  sectional  area  of  the  impacting  object 
and  zero  velocity  everywhere  else,  gives  surprisingly  good 
results. 

3.2  DETAILED  DESCRIPTION  OF  THE  COMPUTATIONAL  PROGRAM 

In  this  section,  a  detailed  description  of  the  computer 
program  used  to  solve  for  the  source  density  distribution,  the 
disturbance  potential,  and  the  disturbance  flow  velocities  along 
with  their  analytical  equations  are  presented. 
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The  Body  Surface  ADproximation 


To  allow  arbitrary  bodies  to  be  considered,  it  is 
required  that  the  body  surface  be  specified  by  a  set  of  points 
in  space  called  input  points.  Then  the  body  surface  is  approxi¬ 
mated  by  a  large  number  of  small  plane  quadrilateral  elements. 

These  elements  are  formed  from  the  original  points  defining  the 
body  surface  as  shown  in  Figure  5 . 

In  the  original  application  of  the  surface  singu¬ 
larity  method  to  FOD  impacts^'  rectangular  surface  elements  were 
used.  However,  in  the  development  of  that  model  no  attempt  was 
made  to  interface  it  with  a  finite  element  structural  analysis 
program,  and  the  impact  surface  was  flat  and  assumed  to  be  rigid. 
Thus,  rectangular  elements  were  preferred  due  to  their  computational 
simplicity.  In  that  model,  the  coordinates  of  the  center  of  each 
element,  and  the  lengths  of  the  long  and  short  sides,  were  read  in 
as  input  data.  The  present  model  is  designed  specifically  to  be 
interfaced  with  finite  element  structural  analysis  computer 


Figure  5.  The  Approximate  Representation  of  the  Body  Surface 


programs  and  is  designed  to  handle  impacts  on  arbitrarily  curved 
surfaces.  Quadrilateral  elements  provide  a  better  approximation 
to  an  arbitrarily  curved  surface  than  rectangular  elements  and 
are  much  easier  to  construct  from  the  finite  element  represen¬ 
tation  of  a  curved  surface.  In  the  model,  it  is  required  that  the 
coordinates  of  all  nodal  points  on  the  impacted  surface  be 
supplied  to  the  loading  model.  Specifically,  the  nodal  points 
in  the  finite-element  program  become  approximate  locations  of 
the  null  point  (defined  to  be  the  point  where  the  quadrilateral 
element  induces  no  velocity  in  its  own  plane)  around  which  the 
plane  quadrilateral  elements  are  constructed  in  the  potential 
flow  analysis.  This  greatly  facilitates  the  coupling  of  the 
finite-difference  loading  model  and  the  finite-element  struc¬ 
tural  code.  Null  points  are  the  points  at  which,  with  the  least 
amount  of  computational  effort  and  with  the  greatest  accuracy, 
the  surface  pressures  can  be  computed  in  the  loading  model.  The 
surface  nodal  locations  are  the  points  where  the  finite  element 
structural  code  requires  the  surface  loads  to  be  defined. 

Input  to  the  loading  model  consists  of  the  coordi¬ 
nates  of  the  set  of  finite  element  nodal  points  defining  the 
surface  of  the  impacted  object.  These  coordinates  are  given  in 
the  reference  coordinate  system,  that  is,  the  coordinate  system 
used  to  describe  the  shape  of  the  impacted  surface  before  impact 
occurs.  Each  nodal  point  is  identified  by  a  pair  of  integers,  m 
ar d  n,  where  n  identifies  the  column  of  points  to  which  it  belongs 
and  m  defines  its  position  in  the  column.  For  a  fan  blade  analysis 
the  undeformed  blade  is  imagined  first  to  be  cut  into  (NS-1)  strips 
by  NS  cuts  (n-lines)  proceeding  from  root  to  tip  with  these  cuts 
orientated  basically  along  chord  l.nes.  Then  the  blade  is  cut 
into  (NC-1)  strips  by  NC  cuts  (m-lines)  proceeding  from  the 
leading  edge  to  trailing  edge.  The  integer  n  thus  ranges  from  1 
to  NS,  and  the  integer  m  ranges  from  1  to  NC.  The  intersections 
of  these  cuts  define  the  initial  position  of  the  surface  nodal 
points  (before  impact) . 
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Quadrilateral  elements  are  computed  from  groups  of 
four  neighboring  nodal  points  which,  for  nonflat  surfaces,  do 
not  generally  lie  in  a  plane.  Consider  a  nodal  point  identified 
by  (m,n)  which  is  not  on  the  first  or  last  m  or  n  lines  (boundary 
lines).  Nodal  points  lying  on  boundary  lines  require  special 
treatment  which  is  discussed  later  in  this  section.  The  first 
step  in  defining  a  quadrilateral  is  to  locate  four  points  around 
the  nodal  point  (m,n)  which  can  be  thought  of  as  "corner"  points 
for  a  surface  with  four  edges  but  which,  in  general,  will  not  be  a 
plane  quadrilateral.  These  four  "corner"  points  are  located 
midway  along  each  of  four  vectors  connecting  the  nodal  point 
(m,n)  with  nodal  points  (m-l,n-l),  (ra+l,n-l),  (m+l,n+l),  and 
(m-l,n+l) .  These  points  are  labeled  1,  2,  3,  and  4  respectively 
as  shown  in  Figure  6  and  will  hereafter  be  called  mid-node 
points.  The  coordinates  of  the  mid-node  points  are: 


Figure  6.  Nodal  Numbering  System  and  Mid-Node  Point  Formulation 
And  Numbering 
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It  should  be  noted  that  once  a  mid-node  point  is 
defined  for  one  element,  it  can  be  used  in  the  formation  of  up 
to  three  neighboring  elements.  Therefore,  it  is  not  necessary 
to  actually  construct  all  four  vectors  shown  in  Figure  6  for 
each  nodal  point. 

Nodal  points  lying  on  boundary  lines  obviously  do 
not  have  four  neighboring  nodal  points.  For  nodal  points  on  the 
leading  edge  of  the  blade  (m*l),  the  node  point  is  not  used  as 
the  null  point.  For  these  points  mid-node  points  1  and  4 
corresponding  to  node  point  (2,n)  are  used  to  define  mid-node 
points  2  and  3  for  the  (l,n)  leading  edge  node  point.  "Corner" 
points  1  and  4  are  defined  by  projecting  half  the  distance  to 
neighboring  leading  edge  node  points  if  they  exist.  If  they  do 
not  exist  [nodes  (1,1)  and  (1,NS)],  then  they  are  artificially 
generated  by  projecting  off  the  blade  along  the  leading  edge  by 
a  distance  equal  to  one-half  the  distance  to  the  leading  edge 
neighbor  which  does  exist. 

Mid-node  points  for  nodes  lying  on  the  other  three 
boundary  lines  are  formed  from  those  mid-node  points  which  they 
share  with  other  nodal  points.  Any  missing  "corner"  points  are 
formed  by  projecting  off  the  surface  of  the  blade  by  an  amount 
which  does  not  exceed  one  half  the  distance  to  its  nearest 
neighboring  nodal  points. 

3.2.2  Formation  of  the  Plane  Quadrilateral  Surface  Element 

The  plane  quadrilateral  surface  elements  are  formed 
from  the  four  appropriate  mid-node  points  as  follows.  First  the 
two  "diagonal"  vectors  T^  and  T2  are  formed  (Figure  7).  Vector 
T^  goes  from  1  to  3  and  vector  T2  goes  from  2  to  4.  In  general, 
these  vectors  are  not  orthogonal  and  do  not  intersect.  Their 
components  are: 


Figure  7.  Formation  of  an  Element  From  Four  Mid-Node  Points 


Tlx  "  X3  “  X1 


Tly  *  *3  *  Y1 


2x 


X4  ”  X2 


L2y 


Y4  “  Y2 


Tl2  “  Z3  “  Z1 
T2x  “  Z4  ~  Z2 


(21) 


The  cross  product,  N,  of  these  vectors  divided  by  its  own  length 
is  taken  as  the  unit  normal  vector,  n,  to  the  plane  of  the  element. 


N 


*2  X  *x 


The  components  of  N  are: 


Nx  “  T2yTlz  “  TlyT2z 
Ny  ’  TlxT2z  '  T2xTlz 
“z  *  T2xTly  -  TlxT2y 

The  components  of  the  unit  normal  vector  n  are: 

n  Nx 
nX  - 

Ny 

ny»-i 

N2 

n2  -  — 
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(22) 


(23) 


(24) 


where  N  -\/n„2  +  **  2  +  N  2  <25) 

To  completely  specify  the  plane  of  the  element  a  point  in  the 
plane  is  also  required.  This  point  is  taken  as  the  point  whose 
coordinates  X/  Y,  2  are  the  averages  of  the  coordinates  of  the 
four  mid-node  points,  i.e., 

X  -  +  x2  +  x3  +  x4) 

7  -  ^(Yx  +  Y2  +  Y3  >►  Y4)  (26) 

7  ■  21  ^  ^  2^) 

Now  the  mid-node  points  will  be  projected  into  the  plane  of  the 
element  along  the  normal  vector.  The  resulting  points  are  the 
corner  points  of  the  plane  quadrilateral  source  element  and 
these,  rather  than  the  mid-node  points,  are  the  points  used  in  all 
calculations.  The  signed  distance  of  the  k-th  raid-node  point 
(k*l,2,3,4)  from  the  plane  is 


dk  »  nx(X  -  Xk)  +  ny { Y  -  Yk)  +  nz(I  -  Zfe)  (27) 

and  the  coordinates  of  the  corner  points  in  the  reference  coor¬ 
dinate  system  are  given  by 


**  +  nx\ 

Yk  +  nydk 

\  *  "A 


(28) 


3.2.3  Formation  of  the  Element  Coordinate  System 


It  is  convenient  to  derive  and  to  use  the  formulas 
for  the  velocities  induced  by  a  quadrilateral  source  element  of 
uniform  strength  at  points  in  space  assuming  the  element  to  lie 
in  a  coordinate  plane.  This  necessitates  constructing  a  coordi¬ 
nate  system  having  two  of  its  axes  in  the  plane  of  the  element. 
Thus,  three  mutually  perpendicular  unit  vectors  are  required,  two 
of  which  are  in  the  plane  of  the  element  and  one  of  which  is 
normal  to  it. 
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The  unit  normal  vector  is  taken  as  one  o£  the  unit 
vectors ,  and  the  unit  vectors  in  the  plane  of  the  element  are 
denoted  by  t^  and  tj.  t^  is  taken  as  divided  by  its  own 
length,  i.e., 


t 


iy 


(29) 


where 

T1  *  V*l*2  *  Tly2  ♦  Tlz2 


(30) 


The  vector  1 2  is  defined  by  tT2  »  n  x  so  that  its  components  are 


t2x  *  ny^lz  "  nztly 
fc2y  "  nzfclx  "  nxtlz 
t2z  “  nxfcly  “  nyfclx 


(31) 


To  transform  the  coordinates  of  points  and  the 
components  of  vectors  between  the  reference  coordinate  system 
and  the  element  coordinate  system,  the  transformation  matrix  is 
required . 


The  elements  of  this  matrix  are  the  components  of 
the  three  basic  unit  vectors,  tT^,  t2,  n. 

The  transformation  matrix  is  thus  the  array 


11  * 

fclx 

a12  "  fcly 

a13 

i 

ft 

K* 

N 

21  " 

t2x 

a22  *  fc2y 

a23 

"  fc2z 

31  " 

nx 

a32  “  ny 

a33 

"  nz 

To  transform  the  coordinates  of  points  from  one 
system  to  the  other,  the  coordinates  of  the  origin  of  the  ele¬ 
ment  coordinate  system  in  the  reference  coordinate  system  are 
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required;  let  these  be  denoted  XQ,  Y0,ZQ.  Then,  if  a  point  has 
coordinates  x',  y' ,  z'  in  the  reference  coordinate  system  and 
coordinates  X,  Y,  Z  in  the  element  coordinate  system,  the 
transformation  from  the  reference  to  the  element  system  is 

X  -  a11(X,-X0)  ♦  a12(Y'-Y0)  ♦  .13(*'-I0) 

Y  "  a2i(X,-V  +  a22(Y'“Yo)  +  a23(Z'-V  <33> 

Z  -  a31(x’-X0)  +  a32(Y’-Yo)  +  a33(z’-Z0) 

while  the  transformation  from  the  element  to  the  reference  system 


X  “  Xo  *  allX  +  a21Y  *  a31Z 
Y  *  Yo  +  a12X  +  a22Y  +  a32Z 
Z'  "  Zo  +  a13X  +  a23Y  +  a33Z 


(34) 


The  origin  is  temporarily  taken  as  the  point  whose  coordinates 
are  the  averages  of  the  four  mid-node  points,  i.e.,  the  point 
with  coordinates  X,  Y,  "Z  in  the  reference  system  ,  and  is  used  to 
find  the  coordinates  of  the  centroid  of  the  area  (Figure  8). 

The  corner  points  are  transformed  into  the  element  coor¬ 
dinate  system  based  on  the  average  point  as  origin.  Their  coor¬ 
dinates  in  the  element  coordinate  system  are  denote  by  nk*'  °* 

Because  they  lie  in  the  plane  of  the  element,  they  have  a  zero S 
coordinate  in  the  element  coordinate  system.  Using  the  above 
transformation  these  coordinates  are; 


£k  *  all^xk”x^  +  a12^Yk“Y^  +  a13^Zk“Z^ 
17k*  *  a21^Xk“X^  *  a22^Yk“Y^  +  a23^Zk“Z^ 


(35) 


The  origin  of  the  element  coordinate  system  is  now  transferred  to 

the  centroid  of  the  area  of  the  quadrilateral.  With  the  average 

point  as  origin,  the  coordinates  of  the  centroid  in  the  element 

4 

coordinate  system  are: 


*o«  "  r)2*)  +£  2*(7?4*  -^1*) 

1  V_  »  -  IV* 


436) 
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Figure  8.  A  Plane  Quadrilateral  Element,  Transfer  of  Origin  from 
Average  Point  to  the  Centroid 


These  are  subtracted  from  the  coordinates  of  the  corner  points  in 
the  element  coordinate  system  based  on  the  average  point  as  ori¬ 
gin  to  obtain  the  coordinates  of  the  corner  points  in  the  element 
coordinate  system  based  on  the  centroid  or  origin,  i.e.. 


Vie'  *  V, 


k  -  1,2, 3, 4 


(37) 


Since  the  centroid  is  to  be  used  as  the  origin  of  the  element 
coordinate  system,  its  coordinates  in  the  reference  coordinate 
system  are  required  for  use  with  the  transformation  matrix. 
These  coordinates  are: 


3.2.4  Determination  of  Null  Point 


It  is  necessary  to  select  a  particular  point  on 
each  quadrilateral  element  where  the  normal  velocity  will  be 
required  to  vanish  and  where  the  flow  velocities  will  be  computed. 
This  point  is  taken  as  the  point  where  the  quadrilateral  induces 
no  velocity  in  its  own  plane.  It  is  designated  the  null  point. 

The  x  and.  y  coordinates  of  this  point,  in  the  element  coordinate 
system,  are  obtained  as  the  solution  of  two  simultaneous,  non¬ 
linear  equations.  These  equations  are 


Vx(X,Y) 
Vy ( X , Y ) 


0 

0 


(39) 


where  V  and  V  are  the  velocity  components  induced  by  a 
x  j 

quadrilateral  source  element  of  unit  source  density.  They  are 
derived  from  the  fundamental  potential  function.  The  general 
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equations  of  the  induced  velocity  components  are: 


*  *3irlos\v??3: 

♦  vx/3"4'*34 


V2wr2+r3'd23 

d23 


r2+r3+d23 


r3+r4+d34 


/r4+rl*d41 
+  -j — -log1  4  1  4i 


‘41 


t: 


r1+r2-d12 


*243,n7r2+r3-d23 

d23  lr2+r3+d23 


♦  £riilog/r3-4-d34 


'34 


r3+r4+d34 


°41 


4^1^41 
r4+rl+d41 


(40) 


(41) 


tan 


+  tan 


-1/  m12e2*h2> 


Zr. 


-l/m23e3"h3N 
l  Zr,  i 


+  tan 


+  tan 


f-^)  -  “"-l( 

l/m23e2~h2\ 

\  Zr2  / 

l/*n34e3“h3\  .._-l/m34e4_"h4\ 

V  Zc3  )  ’  1  zr4  ) 

l/m41e4“h4\  ^  -l/m41e 

r^~)  tan  / 


(42) 
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and 

rk  mJl*  -Ck)2  +  <*  -^k)2  +  z2  *-1*2, 3, 4 

and  * 

«k  -  Z2  +  (X  -$k)2  k«l, 2,3,4 

hk  *  {Y  “\)(X  *^k}  k-1,2,3,4 


(43) 


(44) 


(45) 

(46) 


,  set  equal  to  the  coordinates  of  the  corner 

points,  which  were  obtained  in  the  previous  section,  these 

equations  are  solved  by  means  of  an  iterative  procedure,  which 

utilizes  analytic  expression  for  the  derivatives  of  Vx  and  V^. 

With  the  notation  (  )  ■  3/3  ,  (  )  •  3/3  ,  the  iterative  proce- 

x  a  y  y 

dure  is  as  follows.  Let  Xp  and  Yp  denote  the  p-th  approximation 
to  the  X  and  Y  coordinates  of  the  null  point  and  let  the  notation 
(  denote  the  quantity  in  brackets  evaluated  at  X  »  X  ,  Y  ■  Y  . 

fir  fir 

Once  the  p-th  approximation  has  been  found,  the  (p+l)-th  approxi¬ 
mation  is  obtained  by  solving  the  following  pair  of  linear 
algebraic  equations  for  Xp+^,  Yp+1. 


with  Z  »  0  and 


-  y 


*  V 


-ly  ip) 


(48) 


!(V  )  J(P)(X 


-  V 


[<V  ]  (P)  (V 

y  y  1  p+i 


-  V 


■  -tVy]  (P) 
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The  first  approximation  is  X  ■  Y  ■  0.  The  iterative  procedure  is 
terminated  when  the  induced  velocity  components  at  the  approxi¬ 
mate  null  point  are  both  less  in  absolute  value  than  a  prescribed 
value.  This  value  is  set  at  0.001. 

3.2.5  Formation  of  the  Vector  Matrix  of  Influence 
Coefficients 

The  velocities  induced  by  the  quadrilateral  source 
elements  at  each  other's  null  points  must  be  computed.  This  is 
done  under  the  assumption  that  the  source  density  on  each  element 
is  of  unit  strength.  The  final  result  of  this  calculation  is  the 
complete  set  of  the  velocities  induced  at  each  null  point  by 
every  quadrilateral  element.  This  array  may  be  thought  of  as  a 
"matrix  of  influence  coefficients,"  the  elements  of  which  are 
vectors  in  three-dimensional  space. 

The  basic  caluclation  is  the  computation  of  the 
velocity  components  induced  at  the  null  point  of  the  i-th  element 
by  a  unit  source  density  distribution  on  the  j-th  element. 

The  coordinates  of  the  i-th  null  point  X_  ' ,  Y  ' , 

np  np 

Znp'  are  transformed  into  the  j-th  element  coordinate  system 

obtaining  Xnp,  Ynp,  Znp.  The  transformation  is  accomplished  by 

means  of  Equation  33.  The  velocity  components  are  evaluated  from 

Equations  40,  41,  and  42.  In  these  formulas  X,  Y,  Z  are  replaced 

by  X _ ,  Y _ ,  Z__  and  5.  ,  ru  are  the  coordinates  of  the  corner 

np  np  np  jc  k 

points  of  the  j-th  element. 

In  evaluating  these  velocity  components,  V  and  V 

x  y 

cause  no  trouble.  The  component  V  requires  special  handling  in 

z 

certain  cases.  As  Znp-MJ,  V^O  if  the  null  point  is  approaching  a 
point  in  the  plane  outside  the  boundaries  of  the  quadrilateral. 
However,  2 ir(Sign  Znp)  if  the  null  point  is  approaching  a  point 
within  the  quadrilateral. 

Due  to  round  off  error  which  occurs  while  making 
the  transformation  from  reference  to  element  coordinate  system, 
znp  na¥  have  small  values  with  either  sign  and  the  calculation  of 
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V  will  give  a  value  different  than  what  it  should  be  (i.e.,  0  or 

2ir ) .  To  avoid  this  error,  the  absolute  value  of  2^  is  tested 

np 

before  velocities  are  computed  and,  if  it  is  less  than  some  small 
prescribed  number  (in  this  case  it  is  the  smaller  diagonal/100), 
which  is,  nevertheless,  large  compared  to  the  expected  round-off 
error,  2np  is  set  equal  to  zero  and  Vz  is  set  equal  to  zero  for 
points  outside  the  quadrilateral  or  equal  to  +2*  for  points  inside 
the  quadrilateral.  Another  situation  that  may  cause  trouble  occurs 
when  the  slope  of  a  side  of  the  quadrilateral  is  infinite.  To 
avoid  difficulties  each  of  the  quantities  (C2“51)»  (C3-C2)»  (5 4-53)' 
and  (£^-€4)  are  tested  to  determine  whether  they  are  zero,  and  if 
any  one  of  them  is  zero,  the  two  inverse  tangents  corresponding 
to  that  side  are  set  equal  to  zero.  It  should  be  mentioned  that 
the  inverse  tangents  in  Equation  42  are  evaluated  in  the  normal 
range  -ir/2  to  +ir/2. 

The  induced  velocity  components  V  ,  V  ,  V  are  in 

x  y  z 

the  element  coordinate  system  and  must  be  transformed  to  obtain 

the  components  V  ' ,  V  ' ,  V  '  in  the  reference  coordinate  system. 

x  y  z 

This  is  done  by  using  Equation  33  where  V  .  V„,  V_  replace  X,  Y, 

x  ^  y  ^  z  ^ 

Z,  respectively,  in  these  equations,  and  Vx',  Vy* ,  Vz'  replace 
(X'-X0),  (Y'-Y0),  (Z’-Zq). 

To  obtain  a  set  of  linear  Igebraic  equations  for 
the  unknown  values  of  the  source  density  on  the  elements,  the 
first  step  is  to  calculate  the  normal  velocities  induced  at  each 
null  point  by  the  various  elements,  each  of  which  is  still 
assumed  to  have  a  unit  source  density. 


The  normal  velocity  induced  at  the  null  point  of 
the  i-th  element  by  a  unit  source  density  on  the  j-th  element  is 
obtained  by  taking  the  dot  product  of  v“T  with  the  unit  normal 
vector  of  the  i-th  element  n^.  vT?  is  defined  as  the  vector 
velocity  induced  at  the  null  point  of  the  i-th  element  by  a 
unit  source  density  on  the  j-th  element.  This  induced  normal 
velocity  is  denoted  A^j.  It  is  given  by 


nixV 


+  nijY 


+  nizV 


(49) 
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The  result  is  a  scalar  matrix  whose  elements  are  the  normal 
velocities  induced  at  the  various  null  points  by  the  various 
quadrilateral  elements  with  unit  source  density.  This  matrix 
is  the  coefficient  matrix  of  the  required  set  of  linear  equations, 
since  multiplying  this  by  the  column  matrix  of  the  unknown  values 
of  the  source  density  on  each  element  gives  a  column  matrix  whose 
elements  are  the  true  normal  velocities  induced  at  the  null 
points  by  the  entire  approximate  body  surface.  The  right  hand 
sides  of  the  linear  equations  are  the  negatives  of  the  normal 
components  of  the  onset  flow  at  the  various  null  points. 

3.2.6  Designation  of  the  Onset  Flow 

The  onset  flow  is  used  to  form  a  right  hand  side 
for  use  with  the  coefficient  matrix.  The  onset  flow  is  designated 
by  the  vector  and  must  be  defined  at  each  null  point  i.  This 
vector  will  be  a  spatially  varying  function  for  flows  of  finite 
extent.  It  also  varies  with  time  for  a  deforming  surface. 
Presently,  the  loading  model  uses  the  simple  onset  flow  model 
described  in  paragraph  3.1.2  modified  to  include  the  effect  of 
the  displacement  velocity  of  the  deforming  blade.  The  components 
of  onset  flow  relative  to  the  blade  are  given  by 


Vt  +  Ux 

Vb  +  Uy 
0. 


(50) 


z  z 

for  all  null  points  which  lie  under  the  projected  area  of  the 
impacting  object  on  the  blade.  The  onset  flow  is  specified  as 
zero  for  all  other  null  points.  The  onset  flow  velocity  components, 
given  by  Equation  50,  are  defined  in  the  reference  coordinate 
system.  For  a  blade  analysis,  the  following  reference  coordinate 
system  is  used. 

The  positive  X  direction  is  in  the  direction  of 
rotation  of  the  rotor,  the  positive  y  direction  is  looking  for¬ 
ward  along  the  axis  of  the  engine,  and  the  Z  direction  is  in  the 
plane  of  rotation  of  the  rotor  with  positive  Z  going  from  root 
to  tip.  The  X,  Y,  and  Z  directions  are  shown  in  Figure  2.  In 
Equation  50,  Ux,  U  and  Uz  represent  the  local  deformation 


velocity  components  of  the  blade  and  are  supplied  to  the  loading 
model  by  the  structural  analysis  program. 

More  sophisticated  onset  flow  models  have  been 
developed  but  have  not  as  yet  been  incorporated  into  the  loading 
model.  These  more  sophisticated  onset  flow  models  attempt  to 
take  account  the  spreading  of  the  jet  as  the  impact  surface  is 
approached  and/  thereby/  give  a  better  description  of  the  pressure 
distribution  at  the  edges  of  the  impact  area.  In  Section  4/  a 
more  sophisticated  onset  flow  model  is  presented  which  is  based 
on  two-dimensional  jet  theory. 

3.2.7  The  Linear  Algebraic  Equations  for  the  Values  of 
the  Surface  Source  Density 

3. 2. 7.1  Formulation  of  the  Equations 


Now  the  values  of  the  surface  source 
density  on  the  elements  will  be  obtained  as  the  solution  of  a 
set  of  linear  algebraic  equations.  Recall  that  the  source  den¬ 
sity  is  assumed  constant  on  each  quadrilateral  element.  Thus, 
there  are  N  unknown  values  of  the  source  density,  where  N  is 
the  number  of  elements  formed  from  the  input  point.  The  total 
normal  velocity  is  required  to  vanish  at  the  null  point  of  each 
element  formed  from  the  nodal  points;  therefore,  there  are  N 
equations  for  the  N  unknown  values  of  the  source  density.  The 
total  velocity  induced  at  the  i-th  null  point  by  all  quadrila¬ 
teral  elements  is 


N. 


(51) 


The  normal  component  of  the  onset  flow  at  the  i-th  null  point 
is  the  dot  product  of  the  onset  flow  vector  and  the  unit  normal 
vector  of  the  i-th  element,  i.e.. 


•ni 


ni 


nixV“x  +  niyV<s8y  +  nizV° 


(52) 


The  total  normal  velocity  at  the  i-th  null  point  is  the  sum  of 
Equations  51  and  52.  Thus  the  requirement  that  the  normal  velo¬ 
city  vanish  at  all  null  points  gives  the  following  set  of  linear 
equations  for  the  values  of  the  source  density 
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(53) 


N 


5 

j3*! 


By  solving  this  set  of  equations  a  complete  set  of  source  den¬ 
sities  is  obtained  for  this  onset  flow. 


3. 2. 7. 2  Solution  of  the  Equations 

If  the  impacted  surface  is  perfectly 
flat,  all  of  the  off-diagonal  members  of  A{  .  are  identically 
zero  and  the  solution  of  A^j  is  trivial;  no  matrix  operations 
are  involved.  For  surfaces  which  are  highly  curved,  the  off- 
diagonal  elements  are  generally  small  compared  to  the  diagonal 
elements  and  matrix  iterative  solution  procedures  such  as  the 

4 

Seidel  procedure  must  be  used.  This  procedure  gives  good  con¬ 
vergence  and  requires  less  computer  time  than  direct  elimination 
methods  particularly  when  more  than  200  elements  are  used.  The 
Seidel  procedure  did  not  prove  to  be  a  very  good  method,  however, 
for  blade  analysis.  Difficulties  with  convergence  were  experi¬ 
enced  with  blade  surfaces  which  were  nearly  flat  before  impact 
but  which  became  slightly  deformed  and  concave  during  the  early 
stages  of  impact.  In  this  situation  many  of  the  off-diagonal 
elements  are  of  the  same  order  of  magnitude  as  diagonal  elements 
and  are  not  all  of  the  same  sign.  Also,  because  of  program 
storage  limitations  imposed  on  the  loading  model,  it  was  not 
feasible  to  consider  more  than  80  to  100  elements.  Therefore, 
a  direct  elimination  method  of  solution  is  used  to  solve  the 
equations.  The  method  used  is  Gaussian  elimination  with  pivotal 
condensation.  A  subroutine  called  CSOLVR,  which  was  developed 
by  one  of  the  authors  (Boehman),  is  used  in  the  loading  model 
computer  program.  This  subroutine  was  originally  developed  to 
solve  difficult  linear  systems  such  as  those  encountered  in 
laminar,  compressible  boundary  layer  stability  theory,  and  in 
large  systems  of  chemical  equilibrium  reaction  equations  solved 
by  the  Newton-Raphson  search  procedure. 


3.2.8  Calculation  of  Total  Flow  Velocities  and  Pressures 

Once  the  values  of  the  surface  source  density  have 
been  found,  the  actual  flow  velocities  at  the  null  points  are 
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calculated  by  multiplying  the  elements  o£  the  "matrices  of 
influence  coefficients,"  which  are  the  induced  velocity  components 
that  were  calculated  assuming  a  unit  value  for  all  source  densi¬ 
ties,  by  the  corresponding  true  values  of  the  source  densities, 
and  summing  ail  such  products  that  are  appropriate  for  the  null 
point  in  question.  To  the  results  of  this  summation  must  be 
added  the  proper  components  of  the  onset  flow.  Recall  that  the 
velocity  induced  at  the  null  point  of  the  i-th  element  by  a 
unit  source  density  on  the  j-th  element  is  the  vector  with 
components  X^j,  Y^j,  z^j  corresponding,  respectively,  to  v  ' , 

V  ' ,  V  ’ .  Let  the  total  flow  velocity  at  the  null  point  of  the 

y  z  -► 

i-th  element  be  denoted  by  the  vector  with  components  v^x, 

V^y,  V^2 .  These  components  are  given  by 

Vix  *  J  Xi  j  ffj  +  V®x 

viy  *  £  <^3  +  v«y  <54) 

j-1 

N 

v.  -  V  z.  •  <r  .  +  v  „ 
xz  13  3  ®2 

These  equations  are  evaluated  for  every  null  point.  The  magni¬ 
tude  of  the  velocity  at  each  null  point  is  then  computed  from 

Vi  ■  *  Viy2  +  Vlz2  <55> 

and  finally  the  pressure  at  each  null  point  is  computed.  The 
loading  pressure  at  any  null  point  is  based  on  an  application  of 
Bernoullis  equation  in  the  form 

?i  -  7  Pb(V^2  -  V.2)  (56) 

where  V  ^  in  this  equation  is  defined  as  the  magnitude  of  the 
onset  flow  evaluated  at  the  coordinates  of  the  center  of  impact. 

With  the  simple  onset  flow  model,  negative 
pressures  will  be  obtained  for  flat  or  nearly  flat  surfaces  at 
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null  points  lying  just  outside  of  the  projected  area  of  the 
slice  on  the  blade.3  This  occurs  because  the  tangential  velocity 
induced  by  an  element  in  its  own  plane  approaches  infinity  at 
the  edges  of  the  element.  When  the  surface  is  flat,  or  nearly 
flat,  the  value  of  the  source  density  for  all  null  points  lying 
outside  of  the  immediate  impact  area  is  zero  so  that  the  net 
tang^.tial  velocity  computed  for  these  null  points  is  heavily 
biased  toward  the  tangential  velocity  induced  at  these  null 
points  by  their  newest  neighbor  elements  which  lie  under  the 
projected  area  of  the  slice  on  the  blade.  In  the  current  ver¬ 
sion  of  the  loading  model,  if  a  negative  pressure  is  computed, 
the  negative  value  is  ignored  and  a  zero  pressure  is  assumed 
instead. 

3.2.9  Summary  Description  of  Coupling  Between  Loading 

Model  and  Structural  Dynamic  Analysis 

In  loading  used  for  any  arbitrarily  shaped  impact 
surface,  the  impact  area  is  divided  into  small  flat  elements, 
and  a  uniform  distribution  of  sources  is  assumed  to  cover  each 
area.  At  the  beginning  of  impact,  an  initial  pressure  distri¬ 
bution  is  computed  for  the  undeformed  blade.  During  an  impact 
in  which  local  deformation  takes  place,  the  deformed  shape  of 
the  impact  zone  is  calculated  in  the  dynamic  structural  analy¬ 
sis.  After  significant  deformation  has  occurred,  the  geometry 
of  the  impact  zone  is  provided  to  the  loading  model.  The  loading 
model  is  then  used  to  calculate  a  new  pressure  distribution. 

As  the  structural  analysis  calculation  proceeds,  the  local  shape, 
the  location,  and  the  velocity  of  the  impact  area  are  updated 
and  passed  to  the  loading  model  at  appropriate  time  intervals. 

The  loading  model,  in  turn,  provides  updated  pressure  distribu¬ 
tion  information  for  the  structural  response  computation.  The 
loading  model  is  fully  interactive  with  the  structural  response 
calculation.  The  duration  of  the  impact  is  computed  by  keeping 
track  of  how  much  of  the  slice  has  been  consumed  during  each 
time  interval. 

The  loading  model  is  capable  of  detailed  inter¬ 
action  with  the  structural  response  model  and  of  dealing  with 
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target  translation,  rotation,  and  local  deformation.  The  load- 
response  coupling  modelled  in  this  formulation  is  capable  of 
accurately  predicting  both  overall  target  response  and  local 
deformation. 

The  principal  limitation  of  the  loading  model,  in 
its  present  form,  is  that  it  does  not  include  transient  effects 
of  shock  wave  formation  and  decay.  In  birds,  porosity  appre¬ 
ciably  reduces  the  peak  pressures  without  significantly  affecting 
steady  flow  pressures ,3 ,5  In  addition,  impact  obliquity  reduces 
the  relative  .importance  of  shock  pressures.  Therefore,  it  is 
not  obvious  that  neglect  of  the  shock  aspects  of  bird  impact  on 
blades  is  a  significant  deficiency.  Nevertheless,  efforts  are 
currently  being  made  to  develop  a  simple,  first  order  model 
for  predicting  the  build  up  of  the  peak  pressure  to  the  Hugoniot 
pressure3,5'  and  the  decay  of  pressure  from  the  peak  pressure  to 
the  steady  flow  pressure  by  the  combined  action  of  release  waves 
in  the  soft  body  and  motion  induced  in  the  target  material  by 
the  shock  wave  transmitted  through  the  target  material. 

Preliminary  work  along  these  lines  shows  that  the 
main  effect  of  shocks  during  impacts  on  thin  blades  is  to  impact 
an  initial  deformation  velocity  to  the  blade  material  exposed 
to  the  slice.  Thus,  one  approach  to  handling  shock  effects  may 
be  to  simply  model  the  initial  impact  process  by  imposing  an 
initial  velocity  boundary  condition  to  the  blade  material  in 
the  structural  analysis.  The  appropriate  initial  velocity  to 
be  imported  can  be  computed  from  a  combination  of  Hugoniot 
pressure  variables,  dimensions  of  the  slice,  thickness  of  the 
blade  material,  and  compressibility  properties  of  the  blade 
material . 
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SECTION  4 

FOD  LOADING  MODEL  COMPUTER  PROGRAM 

In  this  section  the  procedures  used  to  translate  the  theo¬ 
retical  aspects  of  the  surface  singularity  technique  into  a 
working  computer  program  are  presented.  The  input  data  and 
program  options  are  described.  Some  details  are  given  on  how 
the  basic  loading  model  is  interfaced  with  finite  element 
structural  analysis  programs.  Some  results  obtained  with  the 
loading  model  in  the  form  of  pressure  distributions  are  presented 
along  with  comparisons  to  experimental  results. 

Some  further  theoretical  developments  associated  with 
improved  onset  flow  models  are  also  presented.  In  particular , 
working  equations  for  an  onset  flow  model  based  on  two-dimensional 
oblique  jet  impacts  are  developed  and  presented. 

The  loading  model  is  set  up  so  that  non-slicing  impacts 
can  be  considered  as  well  as  slicing  impacts.  No  detailed 
experimental  data,  particularly  pressure  distributions,  are 
available  for  slicing  impacts.  Thus,  a  non-slicing  impact 
capability  was  created  for  the  loading  model  to  validate 
the  output  of  the  loading  model  computer  program.  Also, 
non-slicing  capability  was  desired  for  investigating  FOD 
impacts  on  non-rotating  turbine  engine  components  as  well  as 
aircraft  external  surfaces. 

4.1  REFERENCE  COORDINATE  SYSTEM 

While  in  principle  any  Cartesian  coordinate  system  can  be 
used  with  the  surface  singularity  technique  as  developed  in 
Section  III,  it  was  nevertheless  found  to  be  advantageous  to 
work  with  a  specific  coordinate  system  for  handling  slicing 
impacts . 

The  coordinate  system  used  in  this  report  is  one  commonly 
used  in  the  aircraft  engine  industry  and  is  defined  as  follows: 
the  Z  axis  is  taken  in  the  plane  of  rotation  of  the  rotor  with 
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z  »  0  taken  at  the  axis  of  rotation  with  positive  z  pointing 
from  root  to  tip,  the  Y  direction  is  parallel  to  (but  not 
coincident  with)  the  axis  of  rotation  of  the  rotor  and  is 
positive  when  directed  from  the  rear  of  the  engine  toward  the 
front  or  inlet.  The  X  axis  is  in  the  plane  of  rotation  of  the 
rotor  and  is  perpendicular  to  the  z-axis.  The  positive  x  direction 
is  defined  so  as  to  give  a  right-handed  coordinate  system.  For 
a  rotor  viewed  from  the  rear  of  the  engine  looking  forward, 
positive  rotational  speed  is  defined  as  clockwise  rotation. 

Thus,  for  positive  rotational  speed,  the  positive  x-axis  points 
in  the  direction  of  rotation.  In  short  then,  z  is  the  radial 
direction,  y  is  the  forward  direction,  and  x  is  the  tangential 
direction.  The  origin  of  the  coordinate  system,  except  for 
2  *  0,  is  not  fixed  to  a  specific  location.  With  the  coordinate*' 
system  so  defined,  the  program  user  can, in  most  cases,  input 
a  blade  shape  directly  from  design  drawings  generated  in  the 
aircraft  engine  industry. 

4.2  DESCRIPTION  OF  INPUT  DATA  REQUIREMENTS 

The  primary  input  variables  have  already  been  defined  in 
Sections  II  and  III.  They  are  restated  in  this  section  along 
with  the  symbols  used  in  the  loading  model  to  define  these 
quantities . 

4.2.1  Input  Data  for  Slicing  Impacts 

There  are  three  options  for  slicing  impacts.  The 
input  variable  ISLICE  is  used  to  define  the  type  of  object 
being  sliced.  ISLICE  -  1  denotes  a  bird,  ISLICE  -  2  denotes  an 
ice  sphere,  and  ISLICE  *  3  denotes  an  ice  slab.  The  definitions 
of  the  input  variables  for  a  slicing  impact  of  a  bird  are 
given  in  Table  III.  Most  of  the  impact  variables  defined  in 
Table  III  are  used  in  all  three  options.  For  an  ice  sphere 
only  the  definitions  of  Vfa  and  Wfa  are  changed.  For  an  ice 
sphere,  is  taken  as  the  axial  speed  of  the  sphere  (aircraft 
speed)  and  is  the  mass  of  the  ice  sphere.  For  am  ice  slab  the 
user  must  input  for  the  axial  velocity  of  the  slab  relative 
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TABLE  III 

DEFINITION  OF  INPUT  DATA  FOR  SLICING  BIRD  IMPACT. 


Variables 


Symbol  Used  Symbol  Used 

in  Report _ in  Computer  Program 


No.  of  blades  per  stage 


N 


NBL 


Blade  rotational  speed  (rpm)  n 


RPMY 


Axial  speed  of  the  bird  V. 

(user  generally  specifies 
this  as  the  aircraft 
speed) 

Bird  weight 

Radius  on  rotor  at  which  Z. 

impact  occurs 

Time  at  which  impact  begins 

Coordinates  of  blade  lead¬ 
ing  edge  at  Zi 

Coordinates  of  blade  trail¬ 
ing  edge  at  Z^ 

No.  of  chordwise  cuts  on  the  NC 
blade  for  the  purpose  of 
generating  a  grid  system 

No.  of  spanwise  cuts  on  the  NS 

blade  for  the  purpose  of 
generating  a  grid  system 

Mass  density  of  the  bird 


VB 


WB 

RI 

TIM 

(XL, 

YL) 

(XT, 

YT) 

NC 

NS 

RHOB 


to  the  engine  nacelle  (a  number  generally  less  than  the  aircraft 
speed  but  greater  than  zero) .  Also  for  an  ice  slab  the  user  must 
supply  as  input  the  length  of  the  ice  slab  (BL)  and  the  height 
of  the  slab  (DB);  where  height  denotes  how  much  of  the  span  of 
the  blade  is  to  be  exposed  to  the  slab.  The  thickness  of  the 
slice  taken  out  of  the  slab  is  not  an  input  number  but  is  taken 
to  be  the  maximum  possible  slice  width,  h,  computed  in  equation 


The  coordinates  of  the  leading  and  trailing  edge 
points  defined  in  Table  III  are  used  to  define  the  blade  orienta¬ 
tion  angle  5  defined  in  Figure  2. 
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6  -  tan”1  [ (YL-YT) / (XL- XT) ] 

4.2.2  Input  Data  for  Non-Slicing  Impacts 

For  non-slicing  impacts,  a  right  circular  cylinder 
with  L  ■  2D  is  assumed.  A  value  of  one  for  NHL  should  be  inputed 
to  signal  that  non-slicing  impact  is  occurring.  The  location 
of  the  center  of  impact  (XI,  YI,  ZI)  is  then  read  in  as  input 
data  along  with  the  three  velocity  components  of  the  projectile. 

4.3  MAJOR  PROGRAM  VARIABLE  NAMES 

A  list  of  the  major  program  variable  names  is  given  in 
Appendix  A. 

4 . 4  INTERFACING  DETAILS 

In  the  present  version  of  the  loading  model  program  the 
following  decisions  are  assumed  to  be  made  in  the  structural 
analysis  program: 

(1)  the  time  at  which  a  pressure  distribution  is  to  be 
computed  by  the  loading  model 

(2)  the  number  of  chordwise  cuts  NC  and  spanwise  cuts  NS 
to  be  made  on  the  blade  for  the  purpose  of  generating 
a  grid  system  on  the  impact  area 

(i)  the  surface  nodal  locations  where  pressures  are  to 
be  computed 

Whenever  the  structural  analysis  computer  program  decides 
that  an  updated  pressure  distribution  is  required,  the  above 
information  is  supplied  to  the  loading  model  computer  program 
(called  BPRESS) .  The  displacement  velocities  at  the  surface 
nodal  locations  where  pressures  are  to  be  computed  are  also 
supplied  to  BPRESS.  With  this  information  BPRESS  computes  the 
following:  first,  the  corner  points  of  the  plane  quadrilateral 
elements  are  computed  according  to  the  scheme  outlined  in 
Section  3.2.1.  Then  the  element  on  which  the  center  of  impact 
occurs  is  determined.  This  is  done  by  searching  for  the  one 
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surface  element  which  is  intersected  by  the  line  defining  the 
path  of  the  axis  of  symmetry  of  the  impacting  object.  Next 
the  velocity  of  the  impacting  object  relative  to  the  element 
on  which  the  center  of  impact  occurs  is  computed  taking  into 
account  the  displacement  velocity  of  this  element.  This  relative 
velocity  is  used  to  define  the  instantaneous  dynamic  pressure 
of  the  impact  and  is  also  used  to  compute  the  amount  of  impacting 
object  consumed  during  the  time  interval  between  the  previous 
call  to  BPRESS  and  the  current  call  to  BPRESS  (DLC) .  The  amount 
of  object  which  is  not  yet  consumed  (BLR)  and  the  (current  best 
estimate  of  the  time  at  which  the  entire  object  will  have  been 
consumed  (TERM)  are  computed. 

The  vector  matrix  of  influence  coefficients  are  next 
computed  (Section  3.2.5)  followed  by  computation  of  the  onset 
flow.  The  onset  flow  computation  includes  the  effect  of  dis¬ 
placement  velocities.  The  system  of  linear  equations  which 
determine  the  set  of  source  densities  (Section  3.2.7)  is  then 
solved.  The  final  major  computation  in  BPRESS  is  the  computa¬ 
tion  of  flow  velocities  and  pressure  at  the  surface  nodal  points 
which  are  approximately  equivalent  to  the  null  points.  The 
pressures  and  the  estimated  time  of  impact  durations  are  then 
returned  to  the  structural  analysis  executive  routine. 

4.5  PRESSURE  DISTRIBUTIONS  OBTAINED  WITH  THE  LOADING  MODEL 

4.5.1  Oblique  Impacts  on  Rigid  Surfaces 

At  the  present  time,  steady-state  experimental 
pressure  distributions  are  available  only  for  real  and 
substitute  birds  impacting  on  rigid  flat  plates.  Figure  9 
shows  a  comparison  between  computed  and  measured  steady- 
state  pressure  distributions  for  a  45  degree  impact  angle.  The 
distribution  shown  in  Figure  9  is  along  the  major  axis.  Figure 
10  shows  the  same  type  of  comparison  but  for  a  more  oblique 
impact  (25  degrees) .  Figure  11  shows  pressure  distributions 
along  the  minor  axis  for  45  and  25  degree  impacts.  The 
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computed  results  shown  in  these  three  Figures  are  based  on  the 
simple  onset  flow  model  which  was  described  in  Section  3.1.2. 

The  limitations  of  this  simple  onset  flow  model  are  clearly 
evident  from  these  figures  where  it  is  observed  that  in  the  low 
pressure  region,  near  the  edges  of  the  projected  area  of  the 
bird  on  the  plate,  the  computed  pressure  distribution  is  not 
in  agreement  with  the  measured  distribution.  However,  it  can 
be  seen  that  the  agreement  between  theoretical  and  experimental 
results  is  excellent  over  the  impact  region  where  the  pressures 
are  large.  Many  attempts  have  been  made  in  the  course  of  this 
effort  to  develop  better  onset  flow  models  which  lead  to 
improved  pressure  distributions  at  the  edge  of  the  projected 
area  of  the  projectile  on  the  impact  surface.  Except  for  normal 
or  nearly  normal  impacts,  no  reasonably  simple  onset  flow  model 
has  been  discovered  which  yields  better  overall  results  than  the 
very  simplest  model,  that  is,  the  one  used  to  generate  the 
results  shown  in  Figures  9,  10,  and  11. 


Figure  9. 


Normalized  Steady  Flow  Pressure  Distribution  of 
Nominal  1800  g  Real  Bird  (chicken)  Along  Major 
Axis  at  45  Degree  Impact 


NON  *  OtMENStON  AL  PRESSURE  2Py>v 


Figure  10.  Normalized  Steady  Flow  Pressure  Distribution  of 
Nominal  1800  g  Real  Bird  (chicken)  Along  Major 
Axis  at  25  Degree  Impact 


Figure  11.  Comparison  of  Steady  Flow  Pressure  Distributions 
Along  Minor  Axis  for  Ice;  Real  Birds;  and 
Loading  Model 
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4.5.2.  Normal  Impacts  on  Rigid  Surfaces 


Figure  12  shows  the  pressure  distribution  obtained 
with  the  loading  model  using  the  simple  onset  flow  description; 
Experimental  data  points  are  also  shown  in  Figure  12  along 
with  the  theoretical  pressure  distribution  for  a  two-dimensional 
jet.  The  results  shown  in  Figure  12  clearly  show  that  the 
simple  onset  flow  description  is  not  adequate  for  normal  impacts. 
The  total  force  impacted  to  the  impact  surface  corresponding 
to  the  simple  onset  flow  is  only  about  one-half  of  the  theoretical 
value  of  pV^A  where  A  is  the  cross-sectional  area  of  the  pro¬ 
jectile.  The  pressure  distribution  for  a  two-dimensional  jet 
led  to  a  gross  over-estimate  of  the  total  impact  force. 

An  approximate  theoretical  pressure  distribution  for  an  axisym- 
metric  normal  jet  impact  developed  by  Schach  ^  is  also  shown 
in  Figure  12  along  with  Schach 's  measured  pressure  distribution. 

In  the  next  section  an  improved  onset  flow  model 
based  on  a  two-dimensional  jet  theoretical  solution  is  presented 
which  yields  excellent  agreement  with  experimental  pressure 
distributions . 

4.6  AN  IMPROVED  ONSET  FLOW  MODEL  FOR  NORMAL  IMPACTS 

The  inadequacy  of  the  simple  onset  flow  description  for 
normal  impacts  is  due  to  the  fact  that  the  spreading  of  the 
edge  of  the  fluid  jet  as  it  approaches  the  impact  surface  is 
not  taken  into  account.  For  oblique  impacts  this  is  not  a 
serious  shortcoming  since  the  major  portion  of  the  jet  is 
deflected  from  its  oncoming  path  on  that  portion  of  the  impacted 
area  which  lies  within  the  projected  area  of  the  projectile  on 
the  impact  surface.  For  normal  impacts,  however,  much  of  the 
momentum  transfer  occurs  outside  of  the  projected  area  of  the 
projectile  on  the  impact  surface. 

Figure  13  shows  a  comparison  of  the  jet  boundary  free 
streamline  shapes  for  normal  impact  of  two-dimensional  and 
axisymmetric  jets.  This  figure  shows  that  the  turning  of  an 
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axisymmetric  jet  is  much  more  abrupt  than  the  turning  of  a  two- 
dimensional  jet.  From  Figure  13  it  is  observed  that  the  two- 
dimensional  jet  boundary  streamlines  are  not  altogether  dissim¬ 
ilar  from  the  axisymmetric  jet  if  the  Y  axis  of  the  two-dimensional 
jet  is  shifted  downward  by  an  amount  equal  to  the  jet  half-width. 
This  observation  led  to  the  following  question.  Suppose  the  two- 
dimensional  jet  velocity  field  evaluated  at  Y/a  in  the  vicinity 
of  1  was  used  as  an  onset  flow  description.  Is  it  possible  to 
find  some  value  of  Y/a  that  has  a  velocity  field  which  when  used 
as  an  onset  flow  description  would  yield  a  pressure  distribution 
similar  to  the  axisymmetric  pressure  distribution?  Figure  14 
shows  pressure  distribution  that  correspond  to  various  values 
of  Y/a  (denoted  by  Ll  on  Figure  14) .  From  comparison  of  Figure 
14  to  Figure  12,  it  can  be  seen  that  the  pressure  distribution 
for  Y/a  (Ll)  1.25  is  remarkably  close  to  the  measured  axisym¬ 
metric  pressure  distribution.  Thus  if  the  two-dimensional  jet 
velocity  field  evaluated  at  a  non-dimensional  distance  from  the 
impact  surface  of  Y/a  »  1.25  is  used  as  an  onset  flow  descrip¬ 
tion,  a  reasonably  accurate  pressure  distribution  for  axisym¬ 
metric  normal  impacts  on  rigid  flat  plates  can  be  obtained.  In 
the  present  version  of  the  loading  model  computer  program,  this 
improved  onset  flow  description  is  used  instead  of  the  simple 
onset  flow  description  when  normal  impacts  are  considered.  The 
theoretical  Solution  for  the  velocity  field  of  the  two-dimensional 
jet  impacting  on  a  rigid  flat  plate  at  an  incidence  angle  8  is 
presented  in  Appendix  B,  equations  B-8  and  B-9 .  For  a  normal 
impact  (S  *  90  degrees),  it  was  found  that,  at  Y/a  •  1.25,  the 
normal  component  of  velocity  is  approximately  constant  over  the 
jet  width  and  is  equal  to 

V'  -  0.5  58 

and  the  tangential  velocity  distribution  given  by 

U'  -  0.9  (1  -  e  “(x/a) )  59 

Equations  58  and  59  are  used  to  describe  the  onset  flow  velocity 
field  for  normal  impacts  rather  than  the  exact  two-dimensional 


57 


r 


y 


Figure  14.  Comparison  of  Pressure  Distributions  for 
Normal  Impacts 

velocity  solution,  which  must  be  solved  iteratively. 

When  using  the  loading  model  for  normal  impacts,  the  user 
must  define  the  impact  area  to  be  at  least  1.8  times  the 
projectile  radius  in  order  to  obtain  a  pressure  distribution, 
which  gives  a  total  impulse  within  90  percent  of  the  theoretical 
value . 

One  may  question  the  need  of  using  the  loading  model 
computer  program  for  normal  impacts  of  cylindrical  projectiles. 
Why  not  simply  use  the  pressure  distribution  for  axisymmetric 
jets  given  in  Figure  13.  For  normal  impacts  on  rigid  flat 
surfaces  an  eraperical  equation  for  the  pressure  distribution 
would  certainly  be  a  much  simpler  approach.  However,  for  normal 
impacts  on  deformable  surfaces  one  cannot  use  a  rigid  flat 
plate  pressure  distribution  to  describe  the  loading.  The 
present  improved  onset  flow  model  using  equation  59  to  represent 
the  velocity  tangent  to  the  deformed  surface  has  been  found  to 
correctly  describe  the  loading  effects  produced  by  pocketing 
or  cupping  of  the  impact  area  on  deforming  surfaces. 
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4.7  ONSET  FLOW  MODELING  FOR  OBLIQUE  IMPACTS 

A  large  amount  of  effort  has  been  expended  in  attempts  to 
generate  improved  onset  flow  models  for  oblique  models.  Much  of 
this  effort  has  centered  on  using  two-dimensional  jet  theoret¬ 
ical  velocity  field  solutions  in  one  form  or  another  to  model 
the  onset  flow.  While  it  has  been  possible  to  use  two- 
dimensional  jet  theory  to  develop  onset  flow  description,  which 
yield  pressure  distribution  that  agree  with  the  measured  pressure 
distributions  for  oblique  impacts  on  rigid  flat  plates  than  is 
possible  with  the  simple  onset  flow  model,  these  efforts  have 
not  produced  acceptable  results  for  the  following  reason.  The 
programming  and  computer  time  required  to  compute  the  improved 
onset  flow  velocity  field  becomes  very  large,  to  the  point 
where  the  onset  flow  computation  becomes  the  major  computation. 

The  simple  onset  flow  model  yields  acceptable  results  for 
impact  angles  less  than  or  equal  to  45  degrees .  For  normal  or 
nearly  normal  impacts  (90  and  75  degrees)  the  improved  onset  flow 
model  discussed  in  the  previous  section  has  been  found  to  yield 
acceptable  results.  Thus  at  the  present  time  no  relatively 
simple  onset  flow  model  is  available  for  impacts  over  the  45 
to  60  degree  incidence  angle  range.  At  the  present  time  this 
is  not  a  serious  limitation  to  the  scope  of  the  work  of  the 
contractual  effort  since  angles  of  incidence  in  the  45  to  75 
degree  are  not  encountered  in  high  performance  turbine  engines 
(see  Tables  I  and  II  for  APSI  and  F-101  engine  Bird-Blade 
Incidence  angles) . 

4.8  ADDITIONAL  PROGRAMMING  FOR  NON-SLICING  IMPACTS 

As  was  mentioned  in  Section  4.2.2,  when  the  loading  moci.1 
is  used  for  non-slicing  impacts,  the  location  of  the  center  of 
impact  and  velocity  components  of  the  projectile  are  specified 
by  the  user.  The  methods  used  to  treat  this  type  of  impact 
are  identical  to  those  developed  by  the  authors  of  this  report 
for  treating  impacts  on  aircraft  transparencies  and  are  dis¬ 
cussed  in  Appendix  C. 
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APPENDIX  A 
SYMBOL  DEFINITION 

In  this  appendix  the  symbols  used  in  the  program  will  be 
defined  and  related  to  the  equations  derived  or  listed  in 
Section  3. 

N  :  number  of  elements 

XN (I) , YN ( I ) , ZN ( I )  :  X,  Y,  and  Z  coordinates  of  the 

input  points  defining  the  shape 
of  the  impact  surface . 

UX( I) ,UY( I) ,  UZ( I)  :  X,  Y ,  and  Z  components  of  the 

velocity  at  the  null  point 

X(l/I) , X( 2 , 1 ) , X( 3 r I ) , X( 4 , 1)  :  X  coordinates  of  the  corner  points 

forming  the  element 

Y(1,I) /Y(2,I) ,Y(3,I) ,Y{4,I)  :  Y  coordinates  of  the  corner  points 

forming  the  element 

Z  ( 1 , 1 ) / Z ( 2 , 1 )  , Z ( 3 , 1 )  , Z ( 4 , 1 )  :  Z  coordinates  of  the  corner  points 

forming  the  element 

T1X,T1Y,T1Z  :  components  of  the  diagonal  vector 

joining  corners  1  and  3 

T2X,T2Y,T2Z  :  components  of  the  diagonal  vector 

Tj  joining  corners  2  and  4 

XN1,YN1,ZN1  :  components  of  the  cross  product 

vector  N  of  the  two  diagonal 
vectors 
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UNX,UNY,UNZ  :  components  of  the  unit  normal 

vector  n 

XAV,YAV,ZAV  :  coordinates  of  the  average  point 

of  the  four  corner  points  of  the 
element 

D1,D2,D3,D4  :  signed  projection  distance  of  the 

four  input  points  used  to  form 
an  element  into  the  plane  of  the 
element 

XP1,XP2,XP3,XP4  :  X  coordinates  of  the  corner  points 

of  the  quadrilateral  element  in 
the  reference  coordinate  system 

YPl , YP2 , YP3 , YP4  :  Y  coordinates  of  the  corner  points 

of  the  quadrilateral  element  in 
the  reference  coordinate  system 

ZP1,ZP2,ZP3,ZP4  :  Z  coordinates  of  the  corner  points 

of  the  quadrilateral  element  in 
the  reference  coordinate  system 

T1  :  magnitude  of  the  diagonal  vector 

TP1X,TP1Y,TP1Z  :  components  of  the  unit  vector  t^  in 

the  reference  coordinate  system  used 
to  define  the  element  coordinate  system 

TP2X,TP2Y/TP2Z  :  components  of  the  unit  vector  t2  in 

the  reference  coordinate  system  used 
to  define  the  element  coordinate  system 


ZETAP1 , ZETAP2 , ZCTAP3 , ZETAP4  :  ZETA  coordinates  of  the  corner  points 

of  the  quadrilateral  element  in  the 
element  coordinate  system  using  the 
average  point  as  origin 

ETAP1/ETAP2/ETAP3/ETAP4  :  ETA  coordinates  of  the  corner  points 

of  the  quadrilateral  element  in  the 
element  coordinate  system  using  the 
average  point  as  origin 

ZETAORf ETAOR  :  coordinates  of  the  centroid  in  the 

element  coordinate  system  using 
the  average  point  as  origin 

ZETAl f Z ETA2 , ZETA3 , ZETA 4  :  ZETA  coordinates  of  the  corner 

points  of  the  quadrilateral  ele¬ 
ment  in  the  element  coordinate 
system  using  the  centroid  as 
origin 

ETAl / ETA2 , ETA3 , ETA4  :  ETA  coordinates  of  the  corner 

points  of  the  quadrilateral  ele¬ 
ment  in  the  element  coordinate 
system  using  the  centroid  as 
origin 

XO,YO,ZO  :  coordinates  of  the  centroid  of 

the  quadrilateral  element  in  the 
reference  coordinate  system 

XX,YY,ZZ  :  coordinates  of  the  calculated 

null  point  in  the  element  coor¬ 
dinate  system 
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D12,D23,D34,D41  :  length  o£  the  four  sides  of  the 

quadrilateral  element  which  are 
given  in  Equation  43 

R1,R2,R3,R4  :  quantities  defined  by  Equation  45 

U11,U22,U33,U44  :  the  X  components  of  the  velocity 

induced  by  a  side  of  the  quadrila¬ 
teral  element  at  a  null  point 

XU1,YV1,ZW1  :  X,  Y,  and  Z  components  of  the 

velocity  induced  by  a  quadrila¬ 
teral  element  at  a  null  point 

V11,V22,V33,V44  :  the  Y  components  of  the  velocity 

induced  by  a  side  of  the  quadri¬ 
lateral  element  at  a  null  point 

W11,W22,W33,W44  :  Z  components  of  the  velocity 

induced  by  a  side  of  the  quadri¬ 
lateral  element  at  a  null  point 

quantities  used  to  evaluate  the 

partial  derivatives  of  V  and  V 

x  y 

(Equation  48)  with  respect  to  x 
and  y 

VXX , VXY  :  partial  derivatives  of  the  X  com¬ 
ponent  of  the  induced  velocity  Vx 
with  respect  to  X  and  Y 

VYXfVYY  :  partial  derivatives  of  the  Y  com¬ 
ponent  of  the  induced  velocity 
Vy  with  respect  to  X  and  Y 


R1R2X  r  R2R3X / R3R4X , R4R1X  s 
R1R2Y, R2R3Y, R3R4Y, R4R1Y 
0D12 , DD23 , DD34 , DD41 
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DELX,DELY 

XPNP,YPNP,ZPNP 

XNP, YNP,  ZNP 

ZE21 , ZE32 , ZE43 , ZE14 

ET21 , ET32/ET43 , ET14 

El, E2 , E3 , E4 
H1,H2,H3,H4 
M12 ,M23 ,M34 , M41 

UPX, VPY ,WPZ 

AR(X,J) 


change  in  value  of  the  coor¬ 
dinates  of  the  null  point  pro¬ 
duced  by  one  iteration  in  the 
solution  of  the  non-linear 
equations  for  the  null  point 

coordinates  of  the  null  point  of 
the  quadrilateral  element  in  the 
reference  coordinate  system 

coordinates  of  the  i-th  null 
point  in  the  j-th  element  coor¬ 
dinate  system 

:  X  components  of  the  length  of  the 
sides  of  the  quadrilateral 
element 

:  Y  components  of  the  length  of  the 
sides  of  the  quadrilateral 
element 

:  quantities  defined  by  Equations 
44,  46,  and  47 


:  velocity  components  induced  at 
the  null  point  of  the  i-th  element 
by  a  unit  source  density  on  the 
j-th  element  in  the  reference 
coordinate  system 

:  induced  normal  velocity  at  the 
null  point  of  the  i-th  element  by 
a  unit  source  density  on  the  j-th 
element  (Equation  49) 


VINFX , VINFY , VINFZ  s  components  of  a  uniform  onset 

flow  in  the  reference  coordinate 
system 

ON { I )  :  normal  component  of  a  uniform 

onset  flow  at  the  null  point  of 
the  i-th  element  (Equation  52) 

S(I)  :  source  density  on  the  i-th 
quadrilateral  element 

T1,T2/T3  :  components  of  the  velocity 

induced  at  the  null  point  of  the 
i-th  element  by  a  source  density 
S(J)  on  the  j-th  element  in  the 
reference  coordinate  system 

U1,V1,W1  :  components  of  the  total  flow 

velocity  at  the  null  point  of  the 
i-th  element  in  the  reference 
coordinate  system 

VEL  s  magnitude  of  the  total  flow  velo¬ 
city  at  the  null  point  of  the 
i-th  element 

CP  :  pressure  coefficient  at  a  null 
point 

* 

P(I)  j  pressure  magnitude  at  the  null 
point  of  the  i-th  element 
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APPENDIX  3 

TWO-DIMENSIONAL  ONSET  FLOW  MODELING 


To  find  the  exact  source  density  on  the  quadrilateral 
elements,  the  onset  flow  needs  to  be  specified.  A  two- 
dimensional  onset  flow  was  modeled  and  found  adequate  to  use 
with  axisymmetric  jet  impacting  on  flat  surfaces  at  oblique 
angles.  The  onset  flow  is  assumed  to  be  a  two-dimensional  jet 
flow  bounded  by  free  streamlines,  impacting  on  a  flat  surface. 

A  free  streamline  is  a  streamline  which  separates  fluid  in 
motion  from  fluid  at  rest  and  is  a  line  of  constant  speed  and 
pressure.  To  completely  specify  the  onset  flow,  the  boundary 
of  the  jet,  which  is  composed  of  free  streamlines,  along  with  the 
velocity  field  and  the  stagnation  point  need  to  be  derived. 


B-l  DERIVATION  OF  VELOCITY  FIELD  EQUATIONS 


by 


The  velocity  on  the  free  streamline  is  complex  and  denoted 
v  and  written  as: 


v  =  u  -  i  v 


( B— 1 ) 


Milne-Thomson  expressed  2,  where  2  »  x  +  iy,  in  terms  of  v  ,  of 
two  impinging  jets  A^  and  A2,  meeting  and  branching  off  into  two 
other  streams  B^  and  B2  (see  Figure  B-l)  as: 


where  h. ,  h_,  K, ,  and  K,  denote  the  breadths  at  infinity  of  A.,  A,, 

*  *  ^  ^  i  fl  :  y  ^  ^ 

B^,  and  B2;  and  a^  »  U,  a2  »  Ue1  ,  b^  ■  De*  ,  and  b2  ■  Ue^  .  3 

is  the  angle  between  A^  and  B^,  o i  is  the  angle  between  A^  and  A2, 
and  Y  is  the  angle  between  A^  and  B2.  The  expression  for  2  shows 
that  the  motion  is  reversible.  For  our  purposes,  consider  the 
direct  impact  of  two  jets  with  the  same  asymptote  as  shown  in 
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Figure  B-2.  Assume  that  A^  and  Aj  are  two  uniform  streams,  then 
the  inflow  and  outflow  must  balance  to  preserve  continuity. 


h|  +h2  „K,  +  «2 


(B-3) 


From  conservation  of  momentum  in  the  x  and  y  directions  we 
obtain: 


hj+hgcoscc  -  Kj  cos  p  -  Kg  cos  0 

h2  since  -  Kj  sin  £  -  K2sin  0 


Now  it  is  clear  that  a  symmetrical  solution  must  exist, 

K,  =K2  ,  «  -ir  ,  y  =  2 tt-0 

From  Equations  B-4  and  B-3 


Thus 


Cos#  = 


h,-h2 


K,  +  K2 


and 


hl  h2  "  2K 


1. 


( B-5 ) 


Solving  for  h^  and  hj  we  found  that 

h,  =  K|  (I  +  cos£) 
hg  =  Kj  (I  -  COS  0) 

Substituting  these  values  in  the  expression  of  Z  in  Equation  B-2 
we  get: 


>  %  (|  +COS  /9)log  (l '  ~ (*-cos log (l  + 


(B-6) 


^'°9  (' ' Ue  'i)~  e'S*-#109  (' 
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From  reversibility  and  symmetry,  it  is  obvious  that  solving 
Equation  B-6  for  u  and  v  will  give  us  the  onset  flow  needed. 

Non-dimensional izing  the  velocities  u  and  v  by  the  magni¬ 
tude  of  the  impacting  flow  U,  we  get: 


z=k4 


log 


I  -u  +  I  V 


L I  +u‘~  i  v 


+  cos  P  log 


(l-u  +iv')(l+u,-iv') 


-(cos  >9  “  i  sin  >9 

)  log 

|-(u'“  iv')(cos  >9 -  i  sin  p ) 

-  (cos  >9  i  sin  >9 

')  log 

|-(u'-iv,j(cos  p  +  \  sin 

( B-7 ) 


where  u'«  u/U  and  v'  ■  v/U.  Rearranging  Equation  B-7  and 
equating  it  with  Z  *  x  +  iy  gives: 


X  _  2 

Inr 

r*  ^ 

VH-uf+v’2' 

k,/2  w 

Vd+u^+v12 

(B-C) 


+  cos  P  log 


V(l-u)*+v2' \j(\  +u)2  +v'2 


V(l+v'sin/3-u'  cos  /S^+Cusin^+vcos  p)z 

cos  >9 log V (I  "v'sin  p  cos  p?+  (v1  cos  P  -  u'  sin  pf 
(u'  sin  P  +  v'cos  P) 


-  sin  p 


tarT 


-  tan 


r\ 


LO-u'cos/S+v'sin  /3)J 
(v‘  cos  P  -  u'  sin  p) 


(hu'cos  p-v  sin  P) 
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/ 


ta"'&) + cos  13 


tan1 


v' 

l-u‘ 


(B-9) 


/  u*sin  ff+v'cos  B  \  tQrfl  /v'cos  g-usin  0  \ 
\hucos  g+v'sin  £/  \  l"u'cos  g-v'sing/ 


,  sing 
2 


(l+v'sin  g  “  u' cos  £)2+ (using  +  v'cos  g)2 
(l-v1  sin  g  -u' cos  g)2  +  (v'cos  g-  ul  sin  g)2 


Now  the  velocities  u'  and  v'  can  be  calculated,  for  any  point 
within  the  boundary  with  coordinates  x  and  y,  by  solving  simulta¬ 
neously  the  two  non-linear  Equations  (B-8)  and  (B-9).  These 
equations  are  solved  by  means  of  an  iterative  procedure. 

The  iterative  procedure  is  as  follows.  Let  U*  and 
V'p  denote  the  p-th  approximation  of  the  velocity  components  u' 
and  v'  of  any  point  with  known  coordinates  x  and  y.  The  (p+l)-th 
approximation  is  obtained  by  solving  for  u'p+]_  and  v'p+l  the 
pair  of  linear  algebraic  equations 
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The  iterative  procedure  is  terminated  when  the  functions  6  ^  and 
^2  are  both  less  in  absolute  value  than  a  small  prescribed 
value.  The  iterative  procedure  replaces  the  non-linear  Equations 
( B-8 )  and  (B-9)  by  linear  Equations  (B-10),  whose  coefficients 
are  the  derivatives  of  the  non-linear  functions.  In  Equation  (B- 
10)  <5^  and  ^  are  non-linear  functions  of  U'  and  V'  . 


§  (mV)  -  — - _ — _ L,  f(l"u)2+v'^ 

'  ’  k,/2  -  2  log[(|+uf+v.2 


+±cos/3[log 

2  [  (l+u'2+v,2+2v'sin/3-2u,cos  p 


-  log(! +u’2+ v'2-  2  v'  sin  P  “2ii  cos  p) 


(3-11) 


.  r  2u  sin  £-2(u,2+v‘2)cos  p  sin  p 

-  sin  P  ton  - 3 — « - 5 - 5 - ; - 

.  I  ~  (u  +  v;(  sin  73  -  cos*  P )  -  2  u  cos  p 


,  K  y  2  sin  £  l+u2+v2+2v  sin  p- 2u  cos  P 

"■  kj/2  w  2  Ll+u  tv  -2v  sin£-2u  cos/3 


„  t  -1  2uv  -1  2vcos£-2uv 

+  cos  P  tan 1  — tts  -  tan  - # — — - - 

h(u-v2)  |+u2-v2-2u  cos  P 


+  tan 


2v*  1 

I  -  (u^+v2) 


(B-12) 
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The  convergence  is  fairly  rapid  but  care  should  be  taken  while 
evaluating  the  arc  tangents.  They  should  be  evaluated  between  -  it 
and  77 . 


B-2  BOUNDARY  OF  THE  JET 

Consider  the  flow  described  in  Figure  B-2,  let  the  x-axis 
be  along  the  line  of  impact  and  the  origin  be  the  stagnation 
point  of  the  two  flows  when  they  meet  and  branch  into  a  90  degree 
angle.  If  we  regard  the  streamline  y  *  0  as  a  rigid  barrier,  we 
will  get  a  two-dimensional  jet  flow  impacting  on  an  infinite 
plate  at  an  oblique  angle,  which  is  the  case  we  are  studying 
(Figure  B-3)  and  lets  divide  the  flow  in  two  regions,  Region  1 
from  the  centerline  of  the  jet  toward  the  positive  x-axis  and 
Region  2  toward  the  negative  x-axis.  On  a  free  streamline,  the 
complex  velocity 

V  -  0  e-i9  (B‘13) 


If  we  substitute  this  in  the  expression  for  2  in  Equation  (3-2); 
and  equate  the  real  and  imaginary  parts,  we  get  the  coordinates 
(x,y)  of  a  point  on  the  free  streamlines  expressed  in  terms  of  the 
parameter  0.  Substituting  in  Z  we  get: 


z 


+  h2  eia 


^ ~  h2  a  e  +  K  j  £e  k2ye 
log^sin  —  2°-^  -  k(  e1^  log^sin 


(B-14 ) 


-  k2  e'r  log  ^sin  'j 

Substituting  the  values  of  h^,  h^,  ,  Kj,  a and Y ,  calculated 

from  Equations  B-3,  B-4,  and  B-5  into  Equation  B-14  will  give  us 
the  expression  for  Z  of  interest  in  our  problem,  which  is: 
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-  2cos  P 


sin 


In  region  1  9  is  between  -  ■T  and  -(2* -&) ,  let  X  *  -^-9,  then 
0  s  -*T-x  and  x/2  will  be  between  zero  and  r/2  ninus  8/2.  By 
substituting  9  by  its  value  in  Equation  (3-15)  we  get: 

I 


+  log 


log  f  i  cot  -|-j+cos  P  logM^^j+lcos  /3-i  sin/3j 


sin 

-( 

<ir  PS, 

A  "  2  / 

cos 

Pf-i) 

l+l 

(t-!)] 

-  2  cos  /3 


2 


(B- 


-h  log 


-i  sin 


(t  +  t )-(i— f 


£ 


Using  trigonometric  identities  and  equating  2  to  x  +  iy  we  get 

sin  X 


X  *2* 

—  =  —  1 0  sin  p  -  log  (tan  Hr)*  cos  0  log 

kj/2  tt  I  \  2 


cos  X  +  cos  )3 


Y 

k/2 


2_ 

7T 


r  7 r 

2* _  "2"  cos  /3  +  sin  P  log 


‘1 

[l+*on|-)+ 

(i-tan-g 

(l-tan^/21 

l+tan/3/2! 

'l-tan-|)+! 

\{|-tan£/2) 

CVJ 

03.' 

+ 

(B-15) 


i 


16) 


(3-17) 

|  (3-13 
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For  each  value  of  X  between  zero  and  (tt/2  -  8/2)  we  get  the  coor¬ 
dinates  x  and  y  of  a  point  of  the  free  streamline.  A  better  way 
of  writing  the  equation  of  the  streamline  is  to  calculate  X  in 
function  of  y  and  substitute  it  in  the  expression  for  x.  By 
doing  so  we  get: 


X  ♦  ’■ 
2  =tan 


cot  j  tan  h  ~  ^ 


y-  l+cos  p 


sin 


0 


Y 

where  y  = - 

J  M2 


(B-19) 


p 

and  X  =  {  /3  sin  £-log 


Q 

cot  -g"  tanh 


w_h-\  cos  P 
4 


+cos  j9  log 


sin 


r  -i  r  p  r 

2tan'  cot -p- tanh 

*» . » 


0.  ,„r  ir  /y-l+COS^' 


sin  p 

:os 
sin  £ 


COS 


a  ^ 

[2  tan1 


cot 


0 


tanh 


[f( 


y-i+cos  £ 
sin/3 


-pi 


(3-20) 


+  cos£ 


From  Equation  B-  20  we  can  tell  that  this  free  streamline  has  a 
horizontal  asymptote  which  is 

y  ■  1  -  cos  6 

since  7  tends  to  infinity  as  y  approaches  this  value. 

In  region  2  9  is  between  -(2^-®)  and  -Z* ,  let  X  *  (2ir  -  B)  -9, 
then  9  *  -(2t-3)-X  and  X/2  will  be  between  zero  and  6/2. 

Repeating  the  same  procedure  used  for  region  1  we  get: 
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~r  {log 


,  7r  X  P 

i  cot ( *2"  +  *2  2” 


+cos  £  log 


+^cos  /S  -  i  sin  pj  (p  -  i  +  log 


i  sin 


-cos  (•£■  +  T  -  fi)  ■ 


(B-21) 


-  2  cos  £ 


i/3  ,  /.  .  x  \ 

— +log  (isoy) 


Equating  it  to  Z  *  x  +  iy  we  get: 

V/3-tt)  sin  £-log  tan(-|* +  "f"  “§)  +  cos  /3 


2L 

k/2"  7T 


bg-g- 


-Hog 


sin  ^3  -  -  log  sin  cos(4j:  +  -g-  -  jsj) 


_Y_  s  _2_ 
k,/2~  w 


*|-^)+cos  >9)-  sin  /3  log 


sin  y 


cos 


(i-t-sii 


then  X  in  function  of  7  is: 


X  =  Zcof1  loot  P  + 


wy-i-cos  P\ 
sin#  / 


sin  p 


(3-22) 


(B-23) 


(B-24) 
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this  £ree  streamline  also  has  a  horizontal  asymptote: 

y  =  l+  cos  /3 

As  the  parameter  6  approaches  its  limit,  in  either  region  1  or  2 
the  free  streamline  will  approach  an  asymptote.  In  region  1  as  8 
approaches  -{2ir-B),  x  will  approach  (tt-0).  Let  x  -6-6/  and 

substitute  it  in  Equations  B-16,  then  take  the  limit  as  e-H)  to 
get  the  equation  of  the  asymptote.  Substituting  the  value  of  x  in 
the  expression  for  X  and  Y  we  get: 

as  e-*-0 


Eliminate  In  1/e, 


X  _  z 
k,/2  ~  * 


13  sin 


log 


( 


l  +  tan 


k,  sin£ 


iKl-cosyS) 
2  sin  /3 


(B-26) 


Repeating  the  same  procedure  in  region  2,  we  get  the  equation  of 
the  asymptote  in  this  region.  As  8  approaches  —  ( 2 rr—  3 ) ,  x  will 
approach  zero.  Let  x  *  £  and  substitute  it  in  Equations  3-22  and 
B-23. 
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As  e-*  0 


r^2*-  "§■  “f-  (l+cos  j0)  +  sin/3(logsin/3-log-|-) 


k/2 


2. 

k,/2  ^ 


(B-27) 


(/3-*"1sin  /3  -Hog  tan-^  +cos  £  Jiog-g*  +2  logsin  £-log(-|-)jJ 


Eliminate  In  e/2 


y  p  I  /3 

— = j  (>S  -  tt)  sin  p  +  log  tan 


k  /l 


|  .  7T  Y 

+COS  P  \oq~2  +  2logsm  p  +  ^  sfn  ■  ^  * 


(B-28) 


■T  (I+COS0)  ,  . 

2  ~ sirTjs"  ~ °^ 9n ^ 


The  next  step  will  be  to  find  the  x  distance  from  the  stagnation 
point  to  the  centerline  of  the  jet.  To  do  so  we  have  to  find  the 
intersections  of  the  asymptote  with  the  x  axis  and  add  to  it  the 
length  of  the  projected  radius  on  the  x-axis.  (See  Figure  B-3.) 

This  distance  from  the  stagnation  point  to  the  centerline 
of  the  jet  is  denoted  by  C,  and  is  equal  to 

C  -  a  +  b 
or 

C  -  d  -  b 

a  is  found  by  setting  y  equal  to  zero  in  Equation  3-26,  i.e.. 


P  sin  p  -tog^cot  ^J+  cos  P  - 


7T  (l-COS  P) 

2  sin  o  (b-29) 


-  log |^l  + tan  y)2“(l~tan|f 
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and  b  is  expressed  as 


Then  C  becomes 


(B-31) 

(l-tanf)  ) 

The  other  form  of  C,  i.e,  C  »  d  -  b,  is  used  to  check  the 
correctness  of  Equation  B-31. 


C  = 


+  cos  /3 


2im£+  0sin0Hog(cot|) 

V  U-COS0)  (,  a  i 

'T  -H['+'anT) 


b  = 


2  sin  /3 


(B-30) 


Figure  B-3.  The  Two  Regions  of  the  Flow 
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APPENDIX  C 

APPLICATION  OP  LOADING  MODEL  TO  BIRD  IMPACT 
ON  AIRCRAFT  TRANSPARENCIES 
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APPENDIX  C 


APPLICATION  OF  LOADING  MODEL  TO  BIRD  IMPACT 
ON  AIRCRAFT  TRANSPARENCIES 


Bird  impacts  on  aircraft  transparencies  are  different 
from  bird  impacts  on  other  aircraft  surfaces  because 
transparencies  are  not  rigid  structures  under  bird  loading. 

The  transparency  can  significantly  move  and  deform  during  a 
bird  impact.  Therefore,  it  is  necessary  to  consider  the 
windshield  response.  A  transparency  may,  in  general,  respond 
to  impact  in  two  distinctly  different  modes  which  are  termed 
locally  rigid  and  locally  deforming. 

In  the  locally  rigid  case  the  windshield  does  not  signi¬ 
ficantly  deform  in  the  local  area  of  impact.  See  Figure  C-l. 

The  relative  velocity  and  impact  angle  change  during  the 
impact  process,  which  results  in  significant  changes  in  the 
magnitude  and  direction  of  the  force  and  magnitude  of  the 
pressures  excerted  on  the  windshield. 


BEFORE  IMPACT 


DURING  IMPACT 


WINDSHIELD 


Figure  C-l.  Locally  Rigid  Windshield  Response. 
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In  the  locally  deforming  case  the  local  region  of  impact 
undergoes  significant  deformation  including  local  changes  in 
angle  and  shape.  See  Figure  C-2.  The  windshield  forms  a  pocket 
around  the  bird  which  results  in  greatly  increased  local  loading 
and  deformation. 

In  both  modes  the  potential  flow  model  is  best  fitted  to 
sol  for  the  pressure  distribution  on  the  windshield  during 
impact.  Then,  with  the  use  of  a  structural  program,  deformation 
and  velocity  changes,  for  increments  of  tine,  can  be  calculated. 
This  procedure  can  be  repeated  till  the  bird  is  consumed. 

C-l  ADDITIONAL  INFORMATION  NEEDED  FOR  THE  APPLICATION  OF  THIS 

LOADING  MODEL  TO  TRANSPARENCIES 

Initially  the  known  parameters  are  the  components  of  the 
bird  velocity  and  the  impact  point  coordinates.  Since  the  sur¬ 
face  is  going  to  be  subjected  to  rotation  and  deformation  it  is 
important  to  know  at  all  times  on  what  element  the  impact  has 
occured.  This  will  help  calculating  the  relative  velocity  and 
non-dimensional izing  the  velocities  with  respect  to  this  rela¬ 
tive  velocity.  The  direction  cosines  calculated  from  the  bird 
velocity  components  will  give  us  the  direction  of  the  impact. 


Figure  C-2.  Locally  Deforming  Windshield  Response 
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Let  VBX,  VBY,  and  VB2  be  the  components  of  the  bird  velocity  VR. 


where: 


VR  =  Vvbx2+  vby^vbz2 


then:  cos  a  *  VBX/VR 

COS  0  ■  VBY/VR  (C-l) 

cos  y  »  VBZ/VR 

where  a /  8  ,  y  are  the  angles  that  the  velocity  vector  makes, 
respectively,  with  X,  -Y,  and  Z  axis.  The  next  step  is  to  find 
the  element  where  the  impact  occured.  This  is  done  by  trans¬ 
forming  the  components  of  bird  velocity  and  the  coordinates  of 
the  impact  point  XI,  YI,  ZI  into  the  element  coordinate  system. 
This  is  accomplished  by  using  Equation  33.  Then,  try  to  find  the 
intersection  of  the  line,  passing  through  the  impact  point  with 
the  direction  cosines  found  earlier  and  the  plane  (X-Y)  of  the 
element. 

The  equation  of  the  line  is: 


x  -X'  _  y-y*  z-z* 

cos  a  cos  /3  cos  y 


where  X1,  Y* ,  and  Z1  are  the  coordinates  of  the  impact  point  in 
the  element  coordinate  system. 


The  intersection  point  is  found  by  setting  Z  equal  to  zero. 
The  coordinates  of  the  intersection  point  are: 


X 


•  cosq 

■  cosy 


Y 


=  y*  -  r 


cos/3 

cosy 


( C-3 ) 


The  impact  point  is  on  the  element  if  the  sum  of  the  angles, 
formed  by  joining  the  intersection  point  to  the  corner  points  of 
the  quadrilateral  element,  equal  to  360  degrees.  See  Figure  C-3. 
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The  onset  flow  still  needs  to  be  specified.  It  is  set 
equal  to  zero  outside  the  projection  of  the  bird  on  the 
windshield.  Within  the  projection  it  is  equal  to  the  direction 
cosine  minus  the  velocity  of  the  element,  i.e.. 


■  cos  a  -  Uxi 
v»y  *  cos  3  -  Uyi 
V  *  cos  y  -  Uzi 

where  V0BX,  V  ,  and  Vaj2  are  the  components  of  the  velocity  in  the 
x,  y,  and  z  directions  and  Uxi,  Uyi,  and  Uzi  are  the  components 
of. the  non-dimensional  velocity  of  the  element. 

A  testing  procedure  is  set  to  find  what  elements  are  out¬ 
side  the  projection  of  the  bird  on  the  windshield,  i.e.,  what 
elements  have  an  onset  flow  of  zero. 

This  is  accomplished  by  calculating  the  perpendicular 
distance  between  the  null  point  of  the  element  and  the  line  of 
impact.  If  this  distance  is  greater  than  the  radius  of  the  bird, 
the  element  is  outside  the  projection  of  the  bird  and  the  onset 
flow  is  equal  to  zero.  Let  XI,  YI  be  the  intersection  point  of 
the  line  of  impact  and  the  plane  of  the  element,  then  the  perpen¬ 
dicular  distance  D  is: 

r  - 


Vd,2+d,2+d,2' 


where : 


d^  ■  cosy  *  YI 

d2  -  XI  •  cos  Y  ( C-7 ) 

dj  ■  cos  a  ( -YI )  -  cos  g  (-XI) 


Now  we  have  all  the  information  needed  for  the  application  of  the 
loading  model  to  bird  impacts  on  aircraft  transparencies. 
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